Thierry Bernard

Fabrice Vercruyssen

F Grego

Christophe Hausswirth

R Lepers

Jean-Marc Vallier

Jeanick Brisswalter
email: brisswalter@univ-tln.fr

Effect of cycling cadence on subsequent 3 km running performance in well trained triathletes

come

Effect of cycling cadence on subsequent 3 km running performance in well trained triathletes

uring the last decade, numerous studies have investigated the effects of the cycle-run transition on subsequent running adaptation in triathletes. [START_REF] Millet | Physiological and biomechanical adaptations to the cycle to run transition in Olympic triathlon: review and practical recommendations for training[END_REF] Compared with an isolated run, the first few minutes of triathlon running have been reported to induce an increase in oxygen 2-4 fatigue and/or metabolic load induced by a prior cycling event on subsequent running performance. To the best of our knowledge, few studies have examined the effect of cycling task characteristics on subsequent running performance. [START_REF] Garside | Effects of bicycle frame ergonomics on triathlon 10-km running performance[END_REF][START_REF] Hausswirth | Effects of cycling alone or in a sheltered position on subsequent running performance during a triathlon[END_REF][START_REF] Hausswirth | Effect of two drafting modalities in cycling on running performance[END_REF] Hausswirth et al 15 16 indicated that riding in a continuous uptake (V ~O2 ) and heart rate (HR), an alteration in drafting position, compared with the no draft modality, ventilatory efficiency (V ~E), [START_REF] Hue | The influence of prior cycling on biomechanical and cardiorespiratory response profiles during running in triathletes[END_REF] and haemodynamic modifications-that is, changes in muscle blood flow. [START_REF] Kreider | Cardiovascular and thermal response of triathlon performance[END_REF] Moreover, changes in running pattern have been observed after cycling, such as an increase in stride rate 3 6 and modifications in trunk gradient, knee angle in the nonsupport phase, and knee extension during the stance phase. [START_REF] Hausswirth | Relationships between mechanics and energy cost of running at the end of a triathlon and a marathon[END_REF] These changes are generally related to the appearance of leg muscle fatigue characterised by perturbation of electromyographic activity of different muscle groups. [START_REF] Witt | Coordination of leg muscles during cycling and running in triathlon[END_REF] Recently, from a laboratory study, Vercruyssen et al [START_REF] Vercruyssen | Influence of cycling cadences on subsequent running performance in triathlon[END_REF] reported that it is possible for triathletes to improve the adaptation from cycling to running at an intensity corresponding to Olympic distance competition pace (80-85% maximal oxygen uptake (<V>O 2 MAX)). They showed a lower metabolic load during a running session after the adoption of the energetically optimal cadence (73 rpm) calculated from the V ~O2cadence relation [START_REF] Brisswalter | Energetically optimal cadence vs. freely chosen cadence during cycling: effect of exercise duration[END_REF][START_REF] Coast | Linear increase in optimal pedal rate with increased power output in cycle ergometry[END_REF][START_REF] Marsh | The association between cycling experience and preferred and most economical cadences[END_REF][START_REF] Marsh | Effect of cycling experience, aerobic power and power output on preferred and most economical cycling cadences[END_REF] compared with the freely chosen cadence (81 rpm) or the theoretical mechanical optimal cadence (90 rpm). [START_REF] Neptune | A theorical analysis of preferred pedaling rate selection in endurance cycling[END_REF] Furthermore, Lepers et al [START_REF] Lepers | Effect of cycling cadence on contractile and neural properties of knee extensors[END_REF] indicated that, after cycling, neuromuscular factors may be affected by exercise duration or choice of pedalling cadence. They observed, on the one hand, the appearance of neuromuscular fatigue after 30 minutes of cycling at 80% of maximal aerobic power, and, on the other hand, that the use of a low (69 rpm) or high (103 rpm) cycling cadence induced a specific neuromuscular adaptation, assessed by the variation in RMS/M wave ratio interpreted as the central neural input change.

From a short distance triathlon race perspective characterised by high cycling or running intensities, these observations raise a major question about the effect of neuromuscular significantly reduced oxygen uptake during cycling and improved the performance of a 5000 m run in elite triathletes. In addition, Garside and Doran [START_REF] Garside | Effects of bicycle frame ergonomics on triathlon 10-km running performance[END_REF] showed in recreational triathletes an effect of cycle frame ergonomics: when the seattube angle was changed from 73° to 81°, the performance of the subsequent 10 000 m run was improved-that is, there was a reduction in race time.

Therefore, the aim of this study was to examine in outdoor conditions the effects of different pedalling cadences (within the range 60-100 rpm) on the performance of a subsequent 3000 m track run, the latter depending mainly on both metabolic and neuromuscular factors. 17 18

METHODS

Participants

Nine well motivated male triathletes currently competing at the national level participated in the study. They had been training regularly and competing in triathlons for at least four years. For all subjects, triathlon was their primary activity; their mean (SD) times for Olympic distance and sprint distance triathlons were 120 minutes 37 seconds (3.2) and 59 minutes 52 seconds (3.4) respectively. Mean (SD) training distances a week were 9.1 (1.9) km for swimming, 220.5 (57.1) km for cycling, and 51.1 (8.9) km for running. The mean (SD) age of the subjects was 24.9 (4.0) years. Their mean (SD) body weight and height were 70.8 (3.8) kg and 179 (3.9) cm respectively. The subjects were asked to abstain from exhaustive training throughout the experiment. Finally, they were fully informed of the content of the experiment, and written consent was obtained before all testing, according to local ethical committee guidelines.

Maximal cycling test

Subjects first performed a maximal test to determine V ~O2 MAX and ventilatory threshold. This test was carried out on an electromagnetically braked ergocycle (SRM; Jülich, Welldorf, Germany), 19 20 on which the handle bars and racing seat are fully adjustable both vertically and horizontally to reproduce the positions of each subject's bicycle. No incremental running test was performed in this study, as previous investigations indicated similar V ~O2 MAX values whatever the locomotion mode in triathletes who began the triathlon as their first sport. 21 [22 ] This incremental session began with a warm up of 100 W for six minutes, after which the power output was increased by 30 W a minute until volitional exhaustion. During this protocol, V ~O2 , V ~E, respiratory exchange ratio, and HR were continuously recorded every 15 seconds using a telemetric system collecting gas exchanges (Cosmed K4 , Rome, Italy) previously validated by Hausswirth et al. [START_REF] Hausswirth | The cosmed K4 telemetry system as an accurate device for oxygen uptake measurements during exercise[END_REF] V ~O MAX was determined according to criteria described by Howley et al [START_REF] Howley | Criteria for maximal oxygen uptake: review and commentary[END_REF] that is, a plateau in V ~O2 despite an increase in power output, a respiratory exchange ratio value of 1.15, or an HR over 90% of the predicted maximal HR (table 1). The maximal power output reached during this test was the mean value of the last minute. Moreover, the ventilatory threshold was calculated during the cycling test using the criterion of an increase in V ~E/V ~O with no concomitant increase in V ~E/V ~CO . [START_REF] Wasserman | Anaerobic threshold and respiratory gas exchange during exercise[END_REF]

Cycle-run performance sessions

All experiments took place in April on an outdoor track. Outside temperature ranged from 22 to 25°C, and there was no appreciable wind during the experimental period. Each athlete completed in random order three cycle-run sessions (20 minutes of cycling and a 3000 m run) and one isolated run (3000 m). These tests were separated by a 48 hour rest period. Before the cycle-run sessions, subjects performed a 10 minute warm up at 33% of maximal power. [START_REF] Lepers | Effect of cycling cadence on contractile and neural properties of knee extensors[END_REF] During the cycling bout of the cycle-run sessions, subjects had to maintain one of three pedalling cadences corresponding to 60, 80, or 100 rpm. These cycling cadences were representative of the range of cadences selected by triathletes in competition. 15 26 Indeed, it was recently reported that, on a flat road at 40 km/h, cycling cadences could range from 67 rpm with a 53:11 gear ratio to 103 rpm with a 53:17 gear ratio. [START_REF] Lepers | Effect of pedalling rates on physiological response during an endurance cycling exercise[END_REF] However, 60 rpm is close to the range of energetically optimal cadence values, [START_REF] Marsh | Effect of cycling experience, aerobic power and power output on preferred and most economical cycling cadences[END_REF] 80 rpm is near the freely chosen cadence, 6 8 and 100 rpm is close to the cadence used in a drafting situation. 15 16 According to previous studies of the effect of a cycling event on running adaptation, 2 5 the cycling bouts were performed at an intensity above the ventilatory threshold corresponding to 70% of maximal power output (80% V ~O2 MAX) and were representative of a sprint distance simulation. 15 16 The three cycling bouts of the cycle-run sessions were conducted on the SRM system next to the running track. The SRM system allowed athletes to maintain constant power output independent of cycling cadence. In addition, feedback on selected cadence was available to the subjects via a screen placed directly in front of them.

After cycling, the subjects immediately performed the 3000 m run on a 400 m track. The mean (SD) transition time between the cycling and running events (40.4 (8.1) seconds) was the same as that within actual competition. [START_REF] Millet | Physiological and biomechanical adaptations to the cycle to run transition in Olympic triathlon: review and practical recommendations for training[END_REF] During the running bouts, race strategies were free, the only instruction given to the triathlete being to run as fast as possible over the whole 3000 m.

Measurement of physiological variables during the cycle-run sessions

V ~O2 , V ~E, and HR were recorded every 15 seconds with a K4 RQ . The physiological data were analysed during the cycling bouts at the following intervals: 5th-7th minute (5-7), 9th-11th minute (9-11), 13th-15th minute (13-15), 17th-19th minute (17-19), and every 500 m during the 3000 m run (fig 1).

Measurement of biomechanical variables during the cycle-run sessions

Power output and pedalling cadence were continuously recorded during cycling bout. During the run, kinematic data were analysed every 500 m using a 10 m optojump system (MicroGate, Timing and Sport, Bolzano, Italy). From this system, speed, contact, and fly time attained were recorded every 500 m over the whole 3000 m. The stride rate-stride length combination was calculated directly from these values. Thus the act of measuring the kinematic variables had no effect on the subjects' running patterns within each of the above 10 m optical bands.

Blood sampling

Capillary blood samples were collected from ear lobes. Blood lactate was analysed using the Lactate Pro system previously validated by Pyne et al. [START_REF] Pyne | Evaluation of the lactate pro blood lactate analyser[END_REF] Four blood samples were collected: before the cycle-run sessions (at rest), at 10 and 20 minutes during the cycling bouts, and at the end of the 3000 m run.

Statistical analysis

All data are expressed as mean (SD). The stability of the running pattern was described using the coefficient of variation ((SD/mean) 100) for each athlete. [START_REF] Maruyama | Temporal variability in the phase durations during treadmill walking[END_REF] A two way analysis of variance (cadence period time) for repeated measures was performed to analyse the effects of time and cycling cadence using V ~O2 , V ~E, HR, speed velocity, stride variability, speed variability, stride length, and stride rate as dependent variables.

For this analysis, the stride and speed variability (in %) were analysed by an arcsine transformation. A Newmann-Keuls post hoc test was used to determine differences among all cycling cadences and periods during exercise. In all statistical tests, the level of significance was set at p<0.05.

RESULTS

m performances

In this study, the performance of the isolated run was significantly better than the run performed after cycling (583.0 (28.3) and 631.

Running bouts of cycle-run sessions

Table 2 gives mean values for V ~O2 , V ~E, and HR for the running bouts. The statistical analysis indicated a significant interac-running performance. A classical view is that performance in triathlon running depends on the characteristics of the preceding cycling event, such as power output, pedalling 1 29 tion effect (period time + cycling cadence) on V ~O2 during sub-cadence, and metabolic load. Previous investigations have sequent running (p<0.05). V ~O2 values recorded during the run section of the 60 rpm session were significantly higher than during the 80 rpm or the 100 rpm sessions (p<0.05, table 2). These values represent respectively 92.3 (3.0)% (60 rpm run), 85.1 (0.6)% (80 rpm run), and 87.6 (1.2)% (100 rpm run) of cycle V ~O2 MAX, indicating a significantly higher fraction of V ~O2 MAX sustained by subjects during the 60 rpm run session from 1000 to 3000 m than under the other conditions (p<0.05, fig 3). Changes in stride rate within the first 500 m of the 3000 m run were significantly greater during the 80 and 100 rpm run sessions than during the 60 rpm run session (1.52 (0.05), 1.51 (0.05), and 1.48 (0.03) Hz respectively). No significant effect of cycling cadence was found on either stride variability during the run or blood lactate concentration at the end of the cycle-run sessions (table 2).

DISCUSSION

The main observations of this study confirm the negative effect of a cycling event on running performance when compared with an isolated run. However, we observed no effect of the particular choice of cycling cadence on the performance of a subsequent 3000 m run. However, our results highlight an effect of the characteristics of the prior cycling event on metabolic responses and running pattern during the subsequent run.

Cycle-run sessions v isolated run and running performance

shown a systematic improvement in running performance when the metabolic load of the cycling event was reduced either by drafting position [START_REF] Hausswirth | Effects of cycling alone or in a sheltered position on subsequent running performance during a triathlon[END_REF] or racing on a bicycle with a steep seat-tube angle (81°). [START_REF] Garside | Effects of bicycle frame ergonomics on triathlon 10-km running performance[END_REF] Unlike a 3000 m run which is characterised by neuromuscular and anaerobic factors, 17 18 the improvement in running performance in these previous studies was observed over a variety of long distances (5-10 km) where the performance depends mainly on the capacity of the subject to minimise energy expenditure over the whole race. 1 14 15 29 Therefore one explanation for our results is that minimisation of metabolic load through cadence choice during cycling has a significant effect on the running time mainly during events of long duration. Further research is needed into the effect of cadence choice on total performance for running distances close to those of Olympic and Iron man triathlon events.

However, despite the lack of cadence effect on 3000 m race time, our results indicate an effect of cadence choice (60-100 rpm) on the stride pattern or running technique during a 3000 m run. This difference was mainly related to the higher velocity preferred by subjects immediately after cycling at 80 and 100 rpm and to the lower velocity from 1500 to 2500 m after cycling at high cadences. These results may suggest that the use of a low pedalling cadence (close to 60 rpm) reduces variability in running velocity-that is, one of the factors of running technique-during a subsequent run.

For running speeds above 5 m/s (> 18 km/h) and close to maximum values, the change in stride rate is one of the most

To our knowledge only one study has analysed the effect of important factors in increasing running velocity. In our cycling events on subsequent running performance when compared with an isolated run. [START_REF] Hausswirth | Effects of cycling alone or in a sheltered position on subsequent running performance during a triathlon[END_REF] The study showed, during a sprint distance triathlon (0.75 km swim, 20 km bike ride, 5 km run), a significant difference betweena5 km run after cycling (alone and in a sheltered position) and the run performed study, the significant increase in running speed observed during the first 500 m of the 80 and 100 rpm run sessions was associated with a significantly higher stride rate (1.51-1.52 Hz) than in the 60 rpm run session (1.48 Hz). The relation between stride rate and cycling cadence has been reported by 16 without a prior cycling event (isolated run). The cycling event

Hausswirth et al in elite subjects participating in a sprint discaused an increase in mean 5 km race time (1014 seconds) and a decrease in mean running velocity (17.4 km/h) compared with the isolated run (980 seconds and 18.2 km/h).

Our results are in agreement, showing an impairment in running performance after the cycling event whatever the choice of pedalling cadence. There was an increase in mean running time (631 seconds) and a decrease in mean running velocity (17.2 km/h) compared with the performance in the isolated run (583 seconds and 18.5 km/h). Therefore, one finding of our study is that a prior cycling event can affect running performance over the 3 km as well as the 5 km and 10 km distances. 1 29 One hypothesis to explain the alteration in running performance after cycling could be the high metabolic load sustained by subjects at the end of cycling characterised by an increase in blood lactate concentration (4-6 mmol/l) associated with a high V ~O2 MAX (81-83%) and HR max (88-92%). On the other hand, Lepers et al [START_REF] Lepers | Effect of cycling cadence on contractile and neural properties of knee extensors[END_REF] have recently shown in well trained triathletes a reduction in muscular force relating to both central and peripheral factors-that is, changes in M wave and EMG RMS-after 30 minutes of cycling performed at different pedalling cadences (69-103 rpm). We hypothesise that these modifications of neuromuscular factors associated with increasing metabolic load during cycling could increase the development of fatigue just before running, whatever the choice of pedalling cadence.

Cycling cadences and physiological and biomechanical characteristics of running

Our results show no effect of different cycling cadences (60-100 rpm) commonly used by triathletes on subsequent tance triathlon, indicating a significantly higher stride rate after cycling at 102 rpm (1.52 Hz) than after cycling at 85 rpm (1.42 Hz) for the first 500 m of the run.

These observations suggest that immediately after the cycle stage, triathletes spontaneously choose a race strategy directly related to the pedalling cadence, but this effect seems to be transitory, as no significant differences between conditions were reported after the first 500 m of running. This is in agreement with previous studies in which changes in stride pattern and running velocity were found to occur only during the first few minutes of the subsequent run. 1 3 5 6 Furthermore, the fact that triathletes prefer to run at a high pace after cycling at 80 and 100 rpm seems to confirm different anecdotal reports of triathletes. Most triathletes prefer to adopt a high pedalling cadence during the last few minutes of the cycle section of actual competition. Three strategies may be evoked to characterise the choice of cycling cadence: speeding up in the last part of the cycle stage in order to get out quickly on the run (when elite triathletes compete in draft legal events) [START_REF] Millet | Physiological and biomechanical adaptations to the cycle to run transition in Olympic triathlon: review and practical recommendations for training[END_REF] ; reducing power output and spin to minimise the effects of the bike-run transition; maintaining power output while increasing cadence. However, our results show that such a strategy is associated with higher metabolic cost during the cycling stage and greater instability in running pattern, suggesting that it is not physiologically beneficial for the athlete to adopt high pedalling cadences in triathlon competition.

During our study, cycling at 100 rpm was associated with an increase in metabolic cost as classically observed in previous studies for a high cadence such as an increase in V ~O2 , HR, V ~E, [START_REF] Hagan | Effect of pedal rate on cardiorespiratory responses during continuous exercise[END_REF] and blood lactate concentration. [START_REF] Brisswalter | Energetically optimal cadence vs. freely chosen cadence during cycling: effect of exercise duration[END_REF] At the end of the 100 rpm cycling task, mean blood lactate concentration was 7.0 (2.0) mmol/l, suggesting a high contribution of anaerobic metabolism, 8 whereas it was 4.6 (2.1) mmol/l after cycling at 60 rpm. The effect of pedalling rate on physiological adaptation during prolonged cycling has recently been investigated. 8 13 32 Brisswalter et al [START_REF] Brisswalter | Energetically optimal cadence vs. freely chosen cadence during cycling: effect of exercise duration[END_REF] indicated that cycling at a cadence higher than 95 rpm induces a significant increase in V ~O2 , V ~E, and lactate concentration after 30 minutes of exercise in triathletes.

Moreover, our results show an effect of cycling cadence on aerobic contribution during maximal running performance. The subjects were able to sustain a higher fraction of V ~O2 MAX during the 60 rpm run session-that is, 92%-than during the 80 and 100 rpm run sessions-84% and 87% of V ~O2 MAX respectively-(fig 3). These results suggest that the contribution of the anaerobic pathway 17 is more important after the higher cycling rates (80 and 100 rpm) than after the 60 rpm ride and could lead during a prolonged running exercise to earlier appearance of fatigue caused by metabolic acidosis. 33 34 In conclusion, our results confirm the alteration in running performance after a cycling event compared with an isolated run. The principal aim of our investigation was to evaluate the impact of different pedalling rates on subsequent running performance. No significant effect of cycling cadence was found on 3000 m running performance, despite some changes in running strategies, stride rate, and metabolic contributions. We chose a running distance of 3000 m to analyse the possible effect of neuromuscular fatigue-previously reported after a 30 minute cycling exercise at the same intensity [START_REF] Lepers | Effect of cycling cadence on contractile and neural properties of knee extensors[END_REF] -on running performance when neuromuscular and anaerobic factors make important contributions. 17 18 As the effect observed was not significant, the choice of cadence within the usual range does not seem to influence the performance of a middle distance run. One limiting factor of this study may be the choice of a short exercise duration because an effect of metabolic load reduction during the cycling stage on running performance was previously observed for a run longer than 5000 m. For multidisciplinary activities such as triathlon and duathlon, further applied research on the relation between cycling cadence and performance of the subsequent run is required to evaluate the influence of the practical conditions and constraints of actual competition.

Take home message

Compared with an isolated run, completion of a cycling event impairs the performance of a subsequent run independently of the pedalling cadence. However, running strategy, stride rate, and metabolic contribution seem to be improved by the use of a low pedalling cadence (60 rpm). The choice of cycling cadence may have an effect on the running adaptation during a sprint or short distance triathlon. Much research has been conducted on the effects of cycling on physiological variables measured during subsequent running in triathletes. Few authors, however, have examined the effect of variation in cycling task characteristics on either such variables or overall run performance. This study, examining the effect of different pedalling cadences during a cycle at about 80% V ~O2 MAX on performance within a succeeding 3 km run by well trained male triathletes, adds to the published work in this area.

V Vleck

Chair, Medical and Research Committee of the European Triathlon Union and Senior Lecturer, School of Chemical and Life Sciences, University of Greenwich, London, UK Veronica@vleck.fsnet.co.uk

Figure 1

 1 Figure 1Representation of the three cycle-run sessions. TR, Cycle-run transition; BS, blood samples taken; M 1 -M 4 , measurement intervals during cycling at 5-7, 9-11, 13-15, and 17-19 minutes; M 5 -M 10 , measurement intervals during running at 500, 1000, 1500, 2500, and 3000 m; WU, warm up for each condition.

Figure 2 3

 23 Figure 2 Race strategies expressed as the evolution in running Figure 3 Changes in fraction of V ˙ O MAX (FV ˙ O MAX) sustained by

  Ergonomie et performance sportive, UFR STAPS, Université de Toulon-Var, France C Hausswirth, Laboratoire de physiologie et biomécanique, INSEP, Paris, France R Lepers, Groupe analyse du mouvement, UFR STAPS, Université de Bourgogne, France

.

  ................. COMMENTARY ..................

Table 1

 1 Physiological characteristics of the subjects obtained during a maximal cycling testValues are expressed as mean (SD). V ˙ O 2 MAX, maximal oxygen uptake (ml/min/kg); V ˙ EMAX, maximal ventilation (litres/min); HR max , maximal heart rate (beats/min); VT, ventilatory threshold; MAP, maximal power output (W).

	RQ
	2

Table 2

 2 Mean values for power output and speed, oxygen uptake, expiratory flow, heart rate, blood lactate, and running performance obtained during the cycle-run sessions

		Cycle		Cycle		Cycle
	Parameter	(60 rpm)_	Run	(80 rpm)	Run	(100 rpm)	Run
	Power output (W)/speed (km/h)	275.4 (19.4)	17.3 (1.1)	277.1 (18.6)	17.2 (1.20	277.2 (17.2)	17.1 (1.5)
	Oxygen uptake (ml/min/kg)	55.6 (4.6)	62.8 (7.3)*	55.3 (4.0)	57.9 (4.1)	56.5 (4.3)	59.7 (5.6)
	Expiratory flow (litres/min)	94.8 (12.2)	141.9 (15.9)	98.2 (9.2)	140.5 (14.6)	107.2 (13.0)*	140. 5 (21.8)
	Heart rate (beats/min)	163.5 (9.5)	184.2 (4.6)	166.1 (10.4)	185.8 (3.1)	170.7 (4.7)*	182. 6 (5.0)
	Lactataemia (mmol/l)	4.6 (2.1)	9.0 (1.9)	5.1 (2.1)	9.2 (1.2)	7.0 (2.0)*	9.9 (1.8)
	Stride rate (Hz)		1.48 (0.01)		1.49 (0.01)		1.48 (0.02)
	Running performance (s)		625.7 (40.1)		630.0 (44.8)		637.6 (57.9)
	*Significantly different from the other cycle-run sessions, p<0.05.

  1 (47.6) seconds for the isolated run and mean cycle-run sessions respectively). No significant effect of cycling cadence was observed on subsequent 3000 m running performance. Running times were 625.7 (40.1), 630.0 (44.8), and 637.7 (57.9) seconds for the 60, 80, and 100 rpm run sessions respectively (table2). The mean running speed during the first 500 m (fig2) was significantly lower after the 60 rpm ride than after the 80 and 100 rpm cycling bouts(17.

	Cycling bouts of cycle-run sessions
	During the 20 minutes at 60, 80, and 100 rpm cycling bouts,
	average cadences were 61.6 (2.6), 82.7 (4.3) and 98.2 (1.7)
	rpm respectively. Mean HR and V ~E recorded during the 100
	rpm cycling bout were significantly higher than in other
	cycling conditions. Furthermore, blood lactate concentrations
	were significantly higher at the end of the 100 rpm bout than
	after the 60 and 80 rpm cycling bouts (7.0 (2.0), 4.6 (2.1) and
	5.1 (2.1) mmol/l respectively, p<0.05). Conversely, no effect of
	either pedalling rate or exercise duration was found on V ~O2
	(table 2, p>0.05).

5

(1.1), 18.3 (1.1), and 18.3 (1.2) km/h respectively). In addition, the speed variability (from 500 to 2500 m) was significantly lower during the 60 rpm run session than for the other cycle-run conditions (2.18 (1.2)%, 4.12 (2.0)%, and 3.80 (1.8)% for the 60, 80, and 100 rpm run respectively).

2 velocity during the run bouts (60, 80, 100 rpm). *Significantly different from the running velocity during the 60 rpm run session, p<0.05. subjects during the running bouts (60, 80, and 100 rpm). *Significantly different from the initial period, p<0.05; †significantly different from the other conditions, p<0.05.


Student Member, IEEE Julian Freytes

Gilbert Bergna

Member, IEEE Jon Are Suul

Salvatore Are

François D'arco

Frederic Gruson

Hani Colas

Xavier Saad

Guillaud

Salvatore D ' Arco

François Gruson

Frédéric Colas

Member, IEEE Hani Saad

Member, IEEE Xavier Guillaud

Improving Small-Signal Stability of an MMC With CCSC by Control of the Internally Stored Energy

Keywords: HVDC Transmission, Modular Multilevel Converter, State-Space Modeling, Small-Signal Stability Analysis

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The Modular Multilevel Converter (MMC) is currently the most promising topology for High Voltage DC transmission systems (HVDC) [START_REF] Lesnicar | An innovative modular multilevel converter topology suitable for a wide power range[END_REF]. Several advantages of the MMC compared to other topologies may be enumerated, such as lower losses, modularity, scalability and low harmonic content in the output AC voltage [START_REF] Lesnicar | An innovative modular multilevel converter topology suitable for a wide power range[END_REF], [START_REF] Antonopoulos | On dynamics and voltage control of the modular multilevel converter[END_REF]. Nevertheless, the internal dynamics of the MMC topology makes the modeling, control and stability studies of this converter highly challenging [START_REF] Debnath | Operation, control, and applications of the modular multilevel converter: A review[END_REF].

Without control of the internal dynamics, a three-phase MMC will experience large second harmonic currents circulating between the different phases, and potential resonances between the internal equivalent capacitance and the filter inductor [START_REF] Ilves | Steady-state analysis of interaction between harmonic components of arm and line quantities of modular multilevel converters[END_REF]. Thus, a Circulating Current Suppression Controller (CCSC) is commonly applied for eliminating the double frequency circulating currents [START_REF] Tu | Circulating current suppressing controller in modular multilevel converter[END_REF]. However, recent studies have demonstrated that poorly damped oscillations or even instability associated with the DC-side current can occur for MMCs with conventional CCSC-based control [START_REF] Jamshidifar | Small signal dynamic dq model of modular multilevel converter for system studies[END_REF]- [START_REF] Freytes | Small-signal model analysis of droop-controlled modular multilevel converters with circulating current suppressing controller[END_REF].

Due to the multiple frequency components naturally appearing in the internal circulating currents and the arm capacitor voltages of an MMC [START_REF] Harnefors | Dynamic analysis of modular multilevel converters[END_REF], traditional power-system-oriented approaches for state-space modeling, linearization and eigenvalue analysis cannot be directly applied. Thus, in [START_REF] Chaudhuri | Stability analysis of vector-controlled modular multilevel converters in linear time-periodic framework[END_REF], the stability of an MMC was studied by application of time-periodic system theory (Poincaré multipliers), to demonstrate how the double frequency dq current control loops of the CCSC can make the system unstable. Similar conclusions were obtained in [START_REF] Jamshidifar | Small signal dynamic dq model of modular multilevel converter for system studies[END_REF], [START_REF] Li | Stability of a modular multilevel converter based hvdc system considering dc side connection[END_REF], [START_REF] Li | Harmonic instability in mmc-hvdc converters resulting from internal dynamics[END_REF], by eigenvalue analysis within a modeling framework based on dynamic phasors and harmonic superposition, for separately representing the different frequency components of the internal MMC variables. However, recent modeling efforts have lead to the development of state-space models that avoid the approximation of harmonic superposition. Indeed, the MMC model from [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF] is first expressed in a form that separates the variables into groups with only one steady-state oscillation frequency, before the variables are transformed into their associated Synchronously Rotating Reference Frames (SRRFs). The resulting model can be linearized for application of traditional eigenvalue analysis, as demonstrated in [START_REF] Freytes | Small-signal model analysis of droop-controlled modular multilevel converters with circulating current suppressing controller[END_REF] where only the Classical CCSC strategy was considered.

This paper extends the MMC model from [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF] by including a simplified representation of the DC bus voltage dynamics in a Multi-Terminal DC grid (MTDC), represented by an equivalent capacitance and a power source. A DC voltage droop control, as expected for HVDC operation in MTDC grids [START_REF] Beerten | Identification and small-signal analysis of interaction modes in vsc mtdc systems[END_REF], [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF], is also included in the model.

Furthermore, the adapted state-space model is linearized, and the impact of the equivalent DC-side capacitance and the droop gain on the poorly damped oscillation modes that occur with the Classical CCSC are investigated by eigenvalue analysis. Participation factor analysis is applied to identify the source of these oscillations, which is shown to be mainly the interaction between the uncontrolled DC-side current and the voltage or the energy of the internal capacitance of the MMC. This indicates that the stability and control system performance can be improved by introducing closed loop control of the DC-side current and the total stored energy within the MMC. By introducing such Energy-based control, with a structure simplified from [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF], the stability problems are avoided and good performance of the control system can be ensured for the full range of expected operating conditions. The obtained improvement in the small-signal stability is demonstrated by time-domain simulations as well as by eigenvalue analysis.

II. AVERAGED MMC MODELS FOR MATHEMATICAL ANALYSIS

In this section, the Arm Averaged Model (AAM) from [START_REF] Saad | Modular multilevel converter models for electromagnetic transients[END_REF], [START_REF] Saad | Mmc capacitor voltage decoupling and balancing controls[END_REF] is first presented to describe the basic mathematical relations of an MMC. Then, a time-invariant model that can be linearized for small-signal eigenvalue analysis is developed according to the approach presented in [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF].

A. Continuous time-periodic Arm Averaged Model

The Arm Averaged Model (AAM) of the MMC is recalled in Fig. 1. The model presents for each phase j (j = a, b, c), a leg consisting of an upper and a lower arm. Each arm includes an inductance L arm , an equivalent resistance R arm and an aggregated capacitance C arm [START_REF] Saad | Modular multilevel converter models for electromagnetic transients[END_REF].

v G c v G b v G a R f L f R f L f R f L f i ∆ c i ∆ b i ∆ a i U a i U b i U c i L a i L b i L c i dc i dc v U ma v L ma i U ma Carm v U Ca m U a m U b m L a v dc C dc DC Bus i l i l = P l v dc AC source R arm L arm R arm L arm R arm L arm L arm R arm L arm R arm L arm R arm i U mb Carm v U Cb i U mc Carm v U Cc m U c i L mc Carm v L Cc m L c i L mb Carm v L Cb m L b i L ma Carm v L Ca Figure 1
. Schematic of the reference configuration with MMC connected to a DC bus capacitor A simplified model is assumed for the DC bus dynamics consisting of an equivalent capacitor C dc which emulates the capacitance of the DC cables and other converter stations connected to the grid. Also, in parallel with C dc there is a controlled current source i l whose output power is P l as an equivalent model of the power exchanged in the HVDC system.

The modulated voltages v U mj and v L mj , as well as the currents i U mj and i L mj of each arm j are described as follows:

v U mj = m U j v U Cj , v L mj = m L j v L Cj (1) i U mj = m U j i U j , i L mj = m L j i L j (2)
where m U j (m L j ) is the corresponding insertion index, and v U Cj (v L Cj ) is the voltage across the upper (lower) equivalent arm capacitance C arm . The voltage and current of the equivalent capacitor are related by the following equation:

C arm dv U Cj dt = i U mj , C arm dv L Cj dt = i L mj . (3)
For deriving the current dynamics of the AAM, the modulation indexes m ∆ j and m Σ j as well as modulated voltages v ∆ mj and v Σ mj are introduced [START_REF] Saad | Mmc capacitor voltage decoupling and balancing controls[END_REF]:

m ∆ j def = m U j -m L j , m Σ j def = m U j + m L j (4) v ∆ mj def = (-v U mj + v L mj )/2, v Σ mj def = (v U mj + v L mj )/2 (5)
The MMC currents can be expressed as:

i ∆ j def = i U j -i L j , i Σ j def = (i U j + i L j )/2 (6)
where i ∆ j corresponds to the AC grid current, and i Σ j is the common-mode current flowing through the upper and lower arm. The current i Σ j is commonly referred as "circulating current" or "differential current" [START_REF] Saad | Mmc capacitor voltage decoupling and balancing controls[END_REF]. However, the more general term "common-mode current" is preferred in this paper, since i Σ j is calculated as a sum of two currents. The AC grid current dynamics are expressed as:

L ac eq di ∆ j dt = v ∆ mj -v G j -R ac eq i ∆ j ( 7
)
where R ac eq def = (R arm + 2R f )/2 and L ac eq def = (L arm + 2L f )/2. The common-mode arm currents dynamics are given by:

L arm di Σ j dt = v dc 2 -v Σ mj -R arm i Σ j (8)
Finally, the addition and difference of the terms in (3) yields,

2C arm dv Σ Cj dt = m ∆ j i ∆ j 2 + m Σ j i Σ j ( 9
)
2C arm dv ∆ Cj dt = m Σ j i ∆ j 2 + m ∆ j i Σ j ( 10
)
where

v ∆ Cj def = (v U Cj -v L Cj )/2 and v Σ Cj def = (v U Cj + v L Cj )/2.
With the new definitions, the modulated voltages v ∆ mj and v Σ mj can be expressed as follows:

v ∆ mj = - 1 2 m ∆ j v Σ Cj + m Σ j v ∆ Cj (11) v Σ mj = 1 2 m Σ j v Σ Cj + m ∆ j v ∆ Cj (12)
In steady state, "∆" variables are sinusoidal at the fundamental grid frequency ω, while "Σ" variables contain a sinusoidal oscillation at -2ω superimposed to a DC-component [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF]. Thus, the variables can be classified as summarized in Table I.

i ∆ j , v ∆ mj , m ∆ j , v ∆ Cj i Σ j , v Σ mj , m Σ j , v Σ

Cj

In the following section, a non-linear time-invariant model is obtained from the MMC model in "Σ-∆" representation given by ( 7)- [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF].

B. Non-linear time-invariant model using voltage-based formulations in Σ-∆ representation

This section summarizes the time-invariant model of the MMC with voltage-based formulation as proposed in [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF]. To achieve a time-invariant model, it is necessary to refer the MMC variables to their corresponding SRRFs, following the frequency classification shown in Table I. For generic variables x Σ and x ∆ , time-invariant equivalents are obtained with the Park transformation defined in [START_REF] Beerten | Identification and small-signal analysis of interaction modes in vsc mtdc systems[END_REF] as (bold variables means matrix or vectors):

ω ⇒ x ∆ dqz def = x ∆ d x ∆ q x ∆ z = P ω x ∆ a x ∆ b x ∆ c -2ω ⇒ x Σ dqz def = x Σ d x Σ q x Σ z = P -2ω x Σ a x Σ b x Σ c 1 2 1 2 1 2   (13)
The formulation of the MMC variables such that this initial separation of frequency components can be achieved constitutes the basis for the used modelling approach, as illustrated in Fig. 2. As shown in ( 7)-( 12), the Σ and ∆ quantities are not fully decoupled. This results in time-periodic variables in the equations after applying the above transformations. For the Σ variables, time-periodic terms of 6ω are neglected without compromising the accuracy of the model. Conversely, the zero sequences of the vectors in "∆" variables present time-periodic terms of 3ω which has to be taken into account. This component was modeled by means of an auxiliary virtual variable, 90°s

+ω 3ω v ∆ Cabc P ω v Σ Cabc P -2ω v ∆ Cz T 3ω ⊥ 90 • v ∆⊥90 • Cz v Σ Cdqz v ∆ Cdq v ∆ CZ i ∆ abc P ω i Σ abc P -2ω i Σ dqz i ∆ dq m Σ abc P -2ω m Σ dqz m ∆ abc P ω m ∆ z = 0 m ∆ dq -2ω DC MMC Internal State Variables Control variables
hifted from the real one, and by using a transformation T 3ω at +3ω to achieve a model with a defined equilibrium point.

3ω + ⇒ x ∆ Z def =    x ∆ Z d x ∆ Zq    =    cos(3ω) sin(3ω) sin(3ω) -cos(3ω)    T3ω    x ∆ z x ∆90°z   
Using the above definitions, the MMC dynamics in their "Σ-∆" representation can be reformulated as a state-space model where all states reach constant values in steady-sate [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF]. The resulting system is recalled in the following.

1) AC grid currents: Applying the Park transformation at ω to [START_REF] Li | Stability of a modular multilevel converter based hvdc system considering dc side connection[END_REF], the SRRF representation of the AC side current dynamics i ∆ dq are given as:

L ac eq di ∆ dq dt = -v G dq + v ∆ mdq -R ac eq i ∆ dq -J ω L ac eq i ∆ dq ( 14
)
where v G dq is the grid voltage at the point of interconnection and J ω is the cross-coupling matrix at the fundamental frequency, as defined in [START_REF] Saad | Modular multilevel converter models for electromagnetic transients[END_REF],

J ω def = 0 ω -ω 0 . ( 15
)
The AC-side modulated voltage v ∆ mdq is defined in ( 16) as a function of the modulation indexes m ∆ dq and m Σ dqz ,

v ∆ mdq = 1 4 V ∆ m ∆ dq , m Σ dq , m Σ z . ( 16
)
V ∆ is defined as the following 2 × 5 voltage matrix, where all elements are represented in their associated SRRFs:

V ∆ def =   -2v Σ Cz -v Σ Cd ; v Σ Cq ; -v ∆ Cd -v ∆ CZd ; v ∆ Cq + v ∆ CZq ; -2v ∆ Cd v Σ Cq ; v Σ Cd -2v Σ Cz ; v ∆ Cq -v ∆ CZq ; v ∆ Cd -v ∆ CZd ; -2v ∆ Cq   . (17)
2) Common-mode arm currents: Similarly, applying the Park transformation at -2ω to (8), the dynamics of the commonmode arm currents in their time invariant representation i Σ dq and i Σ z are obtained, shown in (18a).

L arm di Σ dq dt = -v Σ mdq -R arm i Σ dq -2J ω L arm i Σ dq (18a) L arm di Σ z dt = -v Σ mz -R arm i Σ z + v dc 2 (18b)
with v dc representing the voltage at the MMC DC terminals.

The modulated voltages driving the currents i Σ dq and i Σ z are v Σ mdq and v Σ mz . These voltages are defined in [START_REF] Freytes | Smallsignal state-space modeling of an hvdc link with modular multilevel converters[END_REF], as a function of the modulation indexes:

v Σ mdqz = 1 4 V Σ m ∆ dq , m Σ dq , m Σ z . (19)
V Σ corresponds to the following 3 × 5 voltage matrix:

V Σ def =      v ∆ Cd + v ∆ CZd -v ∆ Cq + v ∆ CZq 2v Σ Cz 0 2v Σ Cd -v ∆ Cq -v ∆ CZq -v ∆ Cd + v ∆ CZd 0 2v Σ Cz 2v Σ Cq v ∆ Cd v ∆ Cq v Σ Cd v Σ Cq 2v Σ Cz      . (20)
3) Arm capacitor voltages sum: Applying the Park transformation at -2ω to (9), the time invariant dynamics of the voltage sum vector v Σ Cdqz can be expressed by (21).

C arm dv Σ Cdqz dt = i Σ mdqz - -2J ω 0 2×1 0 1×2 0 C arm v Σ Cdqz ( 21
)
with i Σ mdqz representing the modulated current as defined in (22), as a function of the modulation indexes,

i Σ mdqz = 1 8 I Σ m ∆ dq , m Σ dq , m Σ z . ( 22
)
I Σ is the following 3 × 5 time invariant current matrix:

I Σ def =      i ∆ d -i ∆ q 4i Σ z 0 4i Σ d -i ∆ q -i ∆ d 0 4i Σ z 4i Σ q i ∆ d i ∆ q 2i Σ d 2i Σ q 4i Σ z      . ( 23
)
4) Arm capacitor voltages difference: Finally, the steadystate time invariant dynamics of the voltage difference vectors v ∆ Cdq and v ∆ CZ are now recalled. Results are obtained by applying the Park transformation at ω and 3ω to [START_REF] Harnefors | Dynamic analysis of modular multilevel converters[END_REF]. For the sake of compactness, the voltage difference vector is defined

as v ∆ CdqZ def = v ∆ Cd , v ∆ Cq , v ∆ CZ d , v ∆ CZq . C arm dv ∆ CdqZ dt = i ∆ mdqZ - J ω 0 2×2 0 2×2 3J ω C arm v ∆ CdqZ (24) with i ∆
mdqZ representing the modulated current as defined in (25), as a function of the modulation indexes,

i ∆ mdqZ = 1 16 I ∆ m ∆ dq , m Σ dq , m Σ z . ( 25
)
I ∆ is the following 4 × 5 time-invariant current matrix:

I ∆ def =     2i Σ d + 4i Σ z -2i Σ q i ∆ d -i ∆ q 2i ∆ d -2i Σ q -2i Σ d + 4i Σ z -i ∆ q -i ∆ d 2i ∆ q 2i Σ d 2i Σ q i ∆ d i ∆ q 0 -2i Σ q 2i Σ d i ∆ q -i ∆ d 0     . (26)

5) DC bus dynamics:

The DC bus dynamics are modelled by (27), where C dc is the cable model terminal capacitance and P l represents the power injection as seen from the MMC station.

C dc

dv dc dt = P l v dc -3i Σ z ( 27
)
An overview of the model structure corresponding to the MMC and DC bus equations is shown in Fig. 3.

Grid currents v G dq vdc i Σ dqz
Eq. ( 14)

Eq. ( 18)

Calculations Calculations

Eq. ( 16)

Calculations v ∆ mdq v Σ mdqz Eq. (19)
Eq. ( 24)

v ∆ CdqZ Calculations Eq. (21) v Σ Cdqz m Σ dqz m ∆ dq Common-mode currents i Σ dqz i ∆ dq DC bus vdc vdc i Σ z ÷ Pl il
Eq. ( 27) This section presents the control scheme for an MMC as shown in Fig. 4. For the AC-side, the MMC control strategy is based on a classical scheme with two cascaded loops. The outer loop controls the active power P ac following a P ac -v dc droop characteristic with gain k d [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF], and reactive power Q ac . Furthermore, v Σ * md and v Σ * mq are the outputs of the CCSC [START_REF] Tu | Circulating current suppressing controller in modular multilevel converter[END_REF], forcing the circulating currents i Σ dq to zero. The gains of the controllers are calculated to achieve a response time τ Σ i of 5 ms and ζ Σ i of 0.7. The output DC current of the converter is left uncontrolled with this strategy, and it is naturally adjusted to balance the AC and DC power flow.

Inputs States Algebraic equations Differential equations v ∆ CdqZ v Σ Cdqz i ∆ dq v Σ Cdqz v ∆ CdqZ v Σ Cdqz v ∆ CdqZ v Σ mdqz v ∆ mdq
ωL ac eq ωL ac eq v G d v G q i ∆ * d i ∆ * q i ∆ d i ∆ q + - - + v ∆ * md v ∆ * mq + + + - + + PI i ∆ PI i ∆ m ∆ dq m Σ dqz Grid currents control v G d i Σ d + - - + v Σ * md v Σ * mq CCSC i Σ q i Σ * d = 0 i Σ * q = 0 v dc 2 2ωLarm 2ωLarm - - - + ÷ + - v * dc vdc ÷ v G d DC voltage droop control -1 k d + P * ac0 Q * ac P * ac v dc 2 3 2 3 v Σ * mz (28) Equation PI i Σ PI i Σ + -
The MMC insertion indexes are calculated directly from the output of the control loops for the ac-side and CCSC as shown in (28).

m Σ d = 2 v Σ * md v dc , m Σ q = 2 v Σ * mq v dc , m Σ z = 2 v Σ * mz v dc (28a) m ∆ d = -2 v ∆ * md v dc , m ∆ q = -2 v ∆ * mq v dc (28b)
where, for this controller, m Σ z is equal to 1 and v dc is the measured DC voltage. Instead of v dc , a fixed value (i.e. the nominal voltage) can be used, but this choice can reduce the stability region of the MMC [START_REF] Freytes | Small-signal model analysis of droop-controlled modular multilevel converters with circulating current suppressing controller[END_REF]. Calculation of the insertion indexes in this way will be referred to as the "Un-Compensated Modulation" (UCM) [START_REF] Diaz | Small-signal state-space modeling of modular multilevel converters for system stability analysis[END_REF], since there is no compensation for the impact of the oscillations in the equivalent arm capacitor voltages on the modulated output voltages.

IV. SMALL SIGNAL STABILITY ANALYSIS OF AN MMC WITH CLASSICAL CCSC A. Model linearization and time domain validation

The non-linear time-invariant model presented in Section II with the control from Section III are connected as shown in Fig. 5. This interconnected model is represented by a subset of differential equation f as expressed in (29), where x represents the states of the system as in (30) and u the inputs as in (31).

MMC + DC bus

ẋ(t) = f (x(t), u(t)) (29) x = [ξ i ∆ dq ξ i Σ dq Controllers i ∆ dq i Σ dqz v Σ Cdqz v ∆ CdqZ MMC v dc DC bus ] ( 30
) u = [v * dc P * ac0 Q * ac i Σ * d i Σ * q Controllers v G d v G q AC grid ] (31)
The non-linear model from (29) can be linearized around a steady-state operating point by means of the Jacobian linearization method, resulting in a Linear Time-Invariant (LTI) representation as expressed in (32) [START_REF] Wang | Modern power systems analysis[END_REF]. It is recalled that each element A ij and B ij of the matrices A and B are related to the equations f as shown in (33).

∆ ẋ(t) = A(x 0 , u 0 )∆x(t) + B(x 0 , u 0 )∆u(t) (32)
A ij = ∂f i (x, u) ∂x j (x0;u0) ; B ij = ∂f i (x, u) ∂u j (x0;u0) (33)
The equilibrium point, defined by (x 0 , u 0 ), is calculated from the direct resolution of the system equations from (29) by setting ẋ(t) equal to zero. The obtained LTI model is used for evaluating small-signal dynamics and stability by eigenvalue analysis.

To validate the developed small-signal model of the MMC with Classical CCSC, results from a time-domain simulation of two different models will be shown and discussed in the following:

1) EMT: The system from Fig. 1 implemented in EMTP-RV with 400 submodules. The MMC is modeled with the socalled "Model # 2: Equivalent Circuit-Based Model" from [START_REF] Saad | Modular multilevel converter models for electromagnetic transients[END_REF]. It is worth noticing that the modulation indexes are transformed from "Σ-∆" and dq frame to "Upper-Lower" and abc frame. 2) LTI: This model represents the linearized time-invariant model of the interconnected system from Fig. 5, implemented in Matlab/Simulink. The operating point corresponds to the nominal ratings.

The main system parameters are listed in Table II. Starting with a DC power transfer of 1 pu (from DC to AC), a step is applied on P l of -0.1 pu at 0.05 s. The reactive power is controlled to zero during the event. Simulation results are gathered in Fig. 6.

The dynamic response of the DC power P dc is shown in Fig. 6(a) (i.e. measured power in the reference model and the calculated power for the linearized model as P dc = 3i Σ z v dc ). The results of the common-mode currents are shown in Fig. 6(b) (only dq components). The EMT model presents oscillations at 6ω in steady state. These oscillations were neglected during the development of the time-invariant model [START_REF] Bergna | State-space modelling of modular multilevel converters for constant variables in steady-state[END_REF]. As seen in the comparisons from Fig. 6(b), the model captures the average dynamics with reasonable accuracy even if the 6 th harmonic components are ignored (notice the scale). For all other variables, there are negligible differences between the different models.

The step applied on P l produces power imbalance in the DC bus, and the MMC reacts with the droop controller and its internal energy to achieve the new equilibrium point. The DC voltage results are shown in Fig. 6(c). Due to the DC voltage droop control (proportional controller with gain k d ), a steady-state error is obtained after the transient. The internal energy of the MMC participates in the dynamics of the DC voltage regulation by discharging its internal capacitors into the DC bus during the transients, as seen in the voltage v Σ Cz from Fig. 6(d). The behavior of v Σ Cz is similar to the DC bus voltage as expected from the discussion in [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF].

i Σ q -LTI i Σ d -LTI i Σ q -EMT i Σ d -EMT i Σ dq × 10 -3 [pu]
(c) v dc [pu] v Σ Cz -LTI v Σ Cz -EMT v Σ Cz [pu]

B. Stability analysis

Since the linearized model from Fig. 5 has been validated, it can be used for small signal stability analysis. The impact of three main parameters influencing the DC voltage dynamics are evaluated: the DC capacitance, the droop gain k d [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF] and the response time of the current loops.

1) Influence of the DC capacitance: In MTDC systems, the value of the equivalent DC capacitance depends on the number of MMC stations connected to the grid as well as the cable lenghts [START_REF] Freytes | Dynamic impact of mmc controllers on dc voltage droop controlled mtdc grids[END_REF]. This value may vary because some converters could be disconnected or the grid topology reconfigured. For this reason, the MMC should be able to operate under different situations that can result from changes of the DC grid topology and parameters. For evaluating the impact of the DC side capacitance on the small-signal stability, a parametric sweep is performed. The electrostatic constant H dc is defined as,

H dc = 1 2 C dc v 2 dcn P n . ( 34
)
The value of H dc is varied from 40 ms down to 5 ms. This last value represents a small capacitance of the DC bus (24, 4µF << (6 × C arm )), which could represent the DC capacitance of a short cable. The first results consider a power direction from DC to AC side of 1 GW (1 pu) of power transfer. Results are shown in Fig. 7(a). In this case, for the selected values the system remains stable. It is known that the converters dynamics depend on the operating point [START_REF] Freytes | Smallsignal state-space modeling of an hvdc link with modular multilevel converters[END_REF]. The same parametric sweep as the previous example is performed with the opposite power transfer direction (i.e. from AC to DC side). The results are shown in Fig. 7(b), demonstrating that the system becomes unstable when the equivalent DC capacitor decreases. 2) Influence of the droop parameter: In this case, the droop parameter k d is varied from 0.2 pu down to 0.05 pu. The considered power direction is from AC to DC since it is the worst case from the previous section. Results are shown in Fig. 8. When lower values of droop are used, the eigenvalues λ 1,2 shift to the right-hand plane (RHP) resulting in unstable behavior.

C. Identification of unstable eigenvalues

As observed in the previous results (Figs. 7, 8 and 9), the system may become unstable due to the same pair of eigenvalues for all cases (λ 1,2 ). For understanding the origin of these eigenvalues, participation factor analysis is performed for the case from the previous sub-section and the results are shown in Fig. 10, where the considered parameters and operating point are given in Table II. Results from Fig. 10 indicate that the states with the highest participation in the critical mode are i Σ z (i.e. the DC current), v Σ Cz (the state of the MMC which represents the internally stored energy) and v dc (DC voltage). It shows also that the internal circulating currents i Σ dq do not have significant influence on these eigenvalues and neither do the integral part of the controllers (with the chosen bandwidths).

Participation Ratio λ 1,2 0 0.1 0.2 0.3 0.4 v dc v ∆ CZq v ∆ CZ d v ∆ Cq v ∆ Cd v Σ Cz v Σ Cq v Σ Cd i Σ z i Σ q i Σ d i ∆ q i ∆ d ξ i Σ q ξ i Σ d ξ i ∆ q ξ i ∆ d

Controllers

The impact of the proportional gains of the controllers are evaluated by calculating the participation factors for each point from Figs. 9 and the results are shown in Fig. 11. In Fig. 11(a), a similar pattern is observed for the participation factors as in Fig. 10. For fast response times of the CCSC, the dq components of v Σ C participate more on the studied eigenvalues, but the system is unstable as shown in Fig. 9(b). Nevertheless, for realistic values of response times, the most important states are i Σ z , v Σ Cz and v dc , which corresponds to the results in Fig. 10. Since the instabilities identified in the studied cases are due mainly to the uncontrolled output current i Σ z , the natural further step is the explicit control of this current for improving the behavior of the system.

Part. Ratio λ 1,2 τ ∆ i [ms]
v dc v ∆ Cq v ∆ Cd v Σ Cz v Σ Cq v Σ Cd i Σ z (a) Results from participation factor analysis -Eigenvalues λ 1,2 v dc v ∆ Cq v ∆ Cd v Σ Cz v Σ Cq v Σ Cd i Σ z Part. Ratio λ 1,2 τ Σ i [ms] 1.5 2

A. Energy-based controller

The considered control strategy from previous section control two out of three common-mode currents i Σ . The uncontrolled zero-sequence component of i Σ may cause interactions with the DC bus and the internal capacitor voltages, and can potentially make the system unstable. To improve the stability of the studied system, it is proposed to add a DC current control loop (or what is the same, a controller for i Σ z ). In the Classical CCSC strategy from last section, the energy is naturally following the DC bus voltage. The DC current is adjusting itself to achieve an implicit balance of energy and between AC and DC power in steady state. When controlling the DC current, this natural balance is lost so the DC current has to be determined explicitly to regulate the internally stored energy and balance the AC and DC power flow.

1) Inner control loop -Z-sequence Σ current: The design of the controller for i Σ z is based on the second equation of (18a), and a simple PI can be deduced as shown in Fig. 12. For tuning purposes, v Σ mz is supposed to be equal to v Σ * mz .

vdc 2 PI i Σ i Σ z i Σ * z + - v Σ * mz - + Z-seq. Σ Current Controller 1 Larms+Rarm i Σ z vdc 2 + - v Σ mz Equation (18b) ÷ m Σ z Eq. vdc 2 (19) Figure 12. Block diagram of the i Σ z current control
2) Outer control loop -Energy controller: For generating the current reference i Σ * z an outer loop is needed. The proposed strategy is based on the explicit control of the the total stored energy W Σ z on the MMC capacitors C arm given by the power balance between AC and DC sides [START_REF] Saad | Mmc capacitor voltage decoupling and balancing controls[END_REF]. For designing this controller, a model with the explicit relation between the DC current i dc and the total stored energy W Σ z is needed. Assuming P * ac ≈ P ac , a simplified expression of the sum energy dynamics can be defined as [START_REF] Saad | Mmc capacitor voltage decoupling and balancing controls[END_REF]:

dW Σ z dt ≈ P dc -P ac ≈ v dc 3i Σ z i dc -P * ac (35)
The deduced controller structure is shown in Fig. 14. For tuning purposes, the inner i Σ z current controller is considered as a unity gain. Finally, the complete control structure is shown in Fig. 14, where the response time for the total energy τ Σ W is set to 50 ms (i.e. 10 times slower than the inner Σ current loop). The energy reference W Σ * z is set to 1 pu in this paper, for maintaining a constant level of stored energy (corresponding to the rated capacitor voltages). As explained in [START_REF] Freytes | State-space modelling with steady-state time invariant representation of energy based controllers for modular multilevel converters[END_REF], the energy

P * ac PI W Σ W Σ z W Σ * z + - P * dc + Total Energy Controller 1 s W Σ z P * ac - + Equation (35) v dc × i dc 3 i Σ z P dc i Σ * z v dc ÷ 1 3 i * dc i Σ z Ctrl. loop
W Σ z is vdc 2 v dc ÷ 1 3 + - W Σ * z W z + P * ac P * ac i Σ z i Σ * z + - i * dc v Σ * mz - + P * dc Z-seq. Σ Current Controller Total Energy Controller CCSC Grid currents control DC voltage droop control i ∆ * dq v ∆ * mdq i Σ * dq = 0 v Σ * mdq m ∆ dq m Σ dqz v dc PI i Σ PI W Σ Eq. (28)
W Σ z ≈ C arm v Σ Cd 2 2 + v Σ Cq 2 2 + v Σ Cz 2 + ... (36)
...

+ C arm k=d,q,Z d ,Zq v ∆ Ck 2 2

B. Model linearization and time domain validation

In a similar way as in Section IV, the system comprising the non-linear model from Fig. 3 and the controller from Fig. 14 are connected as shown in Fig. 15. The state variables of the new controllers ξ i Σ z and ξ W Σ z , which correspond to the integral part of the PI controllers, are concatenated to the system states x. Also, the energy reference W Σ * z is now included in the inputs vector u. The new interconnected model can be linearized around any operating point for stability analysis. To validate the developed small-signal model of the MMC with Energy based control, results from time domain simulations are shown in Fig. 16. The event and parameters are the same as for Section IV, but it can be observed that the transient behavior of the DC power and voltage are well controlled contrary to the oscillatory behavior from the Classical CCSC strategy (Fig. 6). As shown in Section IV-B, when the Classical CCSC from Fig. 4 was considered, some instabilities were observed with low values of droop gain k d or low equivalent capacitance on the DC side (i.e. low values of H dc ). For demonstrating the stability improvement with the Energy-based controller from Fig. 14, the same parametric sweep is performed as for Fig. 7 and Fig. 8. The results are gathered in Fig. 17. For both situations, it is only shown the case where the power flow is from the DC side to the AC side since it was the case where the instabilities occurred in Section IV-B. However, with this controller, the system presents similar behavior from both power directions.

For the case of the variation of H dc in Fig. 17(a), as well for the variation of k d in Fig. 17(b) it can be clearly observed that stability is guaranteed for the studied cases.

VI. STABILITY COMPARISON

A. DC Capacitance

For comparing the stability improvement, the eigenvalues from Fig. 7(b) and Fig. 17 The stability improvements with the Energy-based controller are highlighted by a time domain simulation. The operating point is the same as for Fig. 18, and results are shown in Fig. 19. Since, for this set of parameters, the configuration of the MMC with Classical CCSC is unstable, the simulation is started with an extra capacitor connected in parallel with C dc for stabilizing the system, which is disconnected at t = 0 s. The frequency of the oscillations corresponds to the frequency of the unstable eigenvalues from Fig. 18.

B. Active power reversal

In this section, an active power reversal from 1 pu to -1 pu is considered. Fixing a DC capacitance with H dc equal to 10ms, and a droop parameter k d equal to 0.1 pu, a parametric sweep for both control strategies is performed and the results are gathered in Fig. 20. For the Classical CCSC in Fig. 20(a) the system is unstable for DC power values lower than -0.15 pu approximately, whereas in Fig. 20(b), i.e. the case with Energy based controller, the system remains stable for all the power range.

VII. CONCLUSIONS

This paper has identified the source of the poorly damped and even possible unstable oscillations that can appear for MMCs relying on a CCSC without control of the internally stored energy. It has been demonstrated by participation factor analysis that the source of these potential stability problems is that the classical CCSC leaves the DC-side current uncontrolled, which leads to potential interaction with the DC bus voltage and the internally stored energy of the the MMC. To avoid these interactions, a control loop for regulating the energy stored in the MMC is introduced. The output of this control loop is the reference for controlling the DC-side current, which is corresponding to the zero-sequence component resulting from the dq-transformation used for the implementation of the classical CCSC. The stability improvement due to the additional control loops is clearly demonstrated by small-signal eigenvalue analysis based on linearization of the presented state-space model, which is accurately representing the internal dynamics of the MMC. The applied state-space modeling approach and the improved dynamic response obtained with the added energyand DC-side current controllers is validated by time-domain simulations in comparison to a detailed EMT model of an MMC with 400 sub-modules per arm. Further developments should consider the dynamic analysis of a more advanced management of the energy in the MMC, i.e. where the energy sum and difference are explicitly controlled.

Figure 2 .

 2 Figure 2. The used modelling approach based on three Park transformations to achieve SSTI state variables

Figure 3 .

 3 Figure 3. MMC and DC bus equations resume

Figure 4 .

 4 Figure 4. MMC droop-controlled with Classical CCSC

Figure 5 .

 5 Figure 3

  P dc[pu]

6 .

 6 Time domain validation Classical CCSC -Step applied on P l of 0.1 pu -EM T : EMTP-RV simulation with detailed converter, LT I: Linear time-invariant state-space model in Simulink

Figure 7 .

 7 Figure 7. Parametric sweep of DC capacitor H dc -DC Operating point v dc0 = 1 pu -k d = 0.1 pu -Classical CCSC

Figure 8 .Figure 9 .

 89 Figure 8. Parametric sweep of Droop parameter k d -DC Operating point v dc0 = 1 pu, P dc0 = -1 pu -Classical CCSC

Figure 10 .

 10 Figure 10. Results from participation factor analysis -Eigenvalues λ 1,2τ ∆ i = 10 ms; τ Σ i = 5 ms; k d = 0.1 pu; H dc = 40 ms

2 Figure 11 .

 211 Figure 11. Participation factors for different response times -Dashed lines correspond to the states that are not listed in the legend -DC Operating point v dc0 = 1 pu, P dc0 = -1 pu

Figure 13 .

 13 Figure 13. Controller design for W Σ z

Figure 14 .

 14 Figure 14. Complete Energy based control -Current and Energy controllers calculated from the dqz components as in (36).

Figure 15 .

 15 Figure 14 controller

Figure 16 .

 16 Figure 16. Time domain validation Energy based controller -Step applied on P l of 0.1 pu -EM T : EMTP-RV simulation with detailed converter, LT I: Linear time-invariant state-space model in Simulink

Figure 17 .

 17 Figure 17. Parametric sweeps for Energy-based controlv dc0 = 1 pu, P dc0 = -1 pu, k d = 0.1 pu -Power flow: DC ⇐ AC

Figure 18 .

 18 Figure 18. Eigenvalues comparison Energy-based and Classical CCSC control -DC Operating point v dc0 = 1 pu, P dc0 = -1 pu, Power flow: DC ⇐ AC k d = 0.1 pu

Figure 19 .

 19 Figure 19. Time domain comparison in EMTP-RV with detailed converters with Classical CCSC and Energy based strategies -Step applied on P l of 0.1 pu -H dc = 14.2 ms

Figure 20 .

 20 Figure 20. Parametric sweep of DC power reversal with Classical CCSC and Energy based strategies -H dc = 10 ms -k d = 0.1 pu Results from Fig. 20 are validated in time domain simulations as shown in Fig. 21. The AC power reference P * ac0 Fig. 4) is also ramped during the transient with the same rate for taking into account the voltage deviation. This represents the case where the power reversal signal is coordinated by a dedicated

  Figure 21. Time domain comparison in EMTP-RV with detailed converters with Classical CCSC and Energy based strategies -H dc = 10 msk d = 0.1 pu

Modular Multilevel Converters connected to Multi-Terminal DC grids in Ecole Centrale de Lille, France. His main research interests include the modeling and control of power electronics for HVDC systems.

Gilbert Bergna received his electrical power engineering degree from the "Simón Bolívar University", in Caracas, Venezuela, in 2008; his "Research Master" from SUPÉLEC (École Supérieure d'Électricité), in Paris, France, in 2010; and a joint PhD degree between SUPÉLEC and the Norwegian University of Science and Technology (NTNU) in 2015. In march 2014 he joined SINTEF Energy Research as a research scientist, and since august 2016 started a post-doctoral fellowship at NTNU working on control and modelling of power electronic systems. Xavier Guillaud (M'04) has been professor in L2EP -Lille since 2002. First, he worked on the modeling and control of power electronic systems. Then, he studied the integration of distributed generation and especially renewable energy in power systems. Nowadays, he is studying the high voltage power electronic converters in transmission system. He is leading the development of an experimental facility composed of actual power electronic converters interacting with virtual grids modelled in real-time simulator. He is involved on several projects about power electronic on the grid within European projects and different projects with French companies. He is member of the Technical Program Committee of Power System Computation Conference (PSCC) and associated editor of Sustainable Energy, Grids and Networks (SEGAN).


Oscar Tellez
email: oscar.tellez@insa-lyon.fr

Samuel Vercraene
email: samuel.vercraene@insa-lyon.fr

Fabien Lehuédé
email: fabien.lehuede@imt-atlantique.fr

Olivier Péton
email: olivier.peton@imt-atlantique.fr

Thibaud Monteiro
email: thibaud.monteiro@insa-lyon.fr

Dial-a-ride problem for disabled people using vehicles with reconfigurable capacity

Keywords: Transportation logistics, optimization, dial-a-ride problem, large neighborhood search meta heuristic, set covering problem

The aim of this paper is to address the dial-a-ride problem with heterogeneous users in which the vehicle capacity can be modified en-route by reconfiguring its internal layout. The work is motivated by the daily transport of children with disabilities performed by a private company based in Lyon Métropole, France. Every day, a fleet of configurable vehicles is available to transport children to medico-social establishments. The objective of this work is then to help route planners with the fleet dimensioning and take reconfiguration opportunities into consideration in the design of routes. Due to the number of passengers and vehicles, real-size instances are intractable for mix-integer programing solvers and exact solution methods. Thus, a large neighborhood search meta-heuristic combined with a set covering component is proposed. The resulting framework is evaluated on real life instances from the transport company.

INTRODUCTION

The standard Dial-a-Ride Problem (DARP) consists in designing vehicle routes in order to serve transportation demands scattered through a geographic area. The global objective is to minimize the transportation cost while satisfying demands.

In contrast to the Pickup and Delivery Problem (PDP), DARP applications concern the transportation of persons. Hence constraints or objectives related to the quality of service should be taken into consideration. In the context of door-to-door transportation of elderly and disabled people, the number of applications has considerably grown recently. Population is aging in developed countries, and many people with disabilities cannot use public transport. As a result, new transport modes, public and private, arise to satisfy their transportation needs. Demand from people with disabilities differ in their need of special equipments such as wheelchair spaces or stretchers thus obliging the use of adapted vehicles. [START_REF] Parragh | Introducing heterogeneous users and vehicles into models and algorithms for the dial-a-ride problem[END_REF] introduced the DARP with heterogeneous users and vehicles, and solved instance with up to 96 user requests. [START_REF] Qu | The heterogeneous pickup and delivery problem with configurable vehicle capacity[END_REF] extended this problem by considering vehicles with configurable capacity. Different categories of users express special needs such as regular seats or wheelchair spaces. These demands are served by configurable vehicles. The goal is to find the most convenient vehicle configuration for each route.

In this paper we present a generalization of the PDP with configurable vehicle capacity. Contrary to the work by [START_REF] Qu | The heterogeneous pickup and delivery problem with configurable vehicle capacity[END_REF], we allow vehicles to be reconfigured enroute on. The other difference is that we determine the fleet dimension instead of starting with limited fleet. We call this variant dial-a-ride problem with reconfigurable vehicle capacity (DARP-RC). Note it is not an heterogeneous variant because only one vehicle type is considered.

The use of hybrid methods or matheuristics has become quite popular in recent years for routing problems. Our solution method has been inspired by the framework of [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with crossdocking[END_REF] to solve the vehicle routing problem with cross-docking. We combine a Large Neighborhood Search (LNS) metaheuristic with a Set Covering Problem (SCP).

The contribution of this paper is therefore to introduce and solve the DARP-RC. Moreover, we compare the results of the combined LNS-SCP approach with pure LNS and adaptive large neighborhood search (ALNS).

INDUSTRIAL CONTEXT

This work is motivated by the daily transport of people with disabilities at Lyon Métropole. One segment of this service is operated by the GIHP company on regular basis. Every day, a fleet of configurable vehicles transport around 500 children from and to Medico-Social Establishments (MSE) for rehabilitative treatment. GIHP serves around 60 MSEs with around 180 adapted vehicles. One of the particularities of the vehicle fleet is the capacity to reconfigure its internal layout to trade off seats by wheelchair spaces as per convenience. 1 For MSEs, transportation is often considered the second biggest expense after wages. As a consequence, optimizing the transport becomes a priority. Every year, the company makes strategic choices in the definition and constitution of this fleet. Then, routing decisions are re-revalued daily by route planners.

These decisions are often taken without help of decision making tools like a vehicle routing software. This is why route planners conceive suboptimal assumptions such as designing separate routes for each MSE or ignore vehicle reconfiguration possibilities. As such, this measures can reduce pooling gains and increase operating cost.

PROBLEM DEFINITION

In the classic DARP, a homogeneous vehicle fleet is assumed. All vehicles have the same single capacity type and are located at a single depot [START_REF] Cordeau | A tabu search heuristic for the static multi-vehicle dial-a-ride problem[END_REF]. The proposed DARP-RC constitutes an extension for DARP problem considering more realistic assumptions such us heterogeneous users (e.g. seats, wheelchairs, stretchers) and vehicles with reconfigurable capacity.

Reconfigurable vehicles

Vehicles with configurable capacity were introduced in [START_REF] Qu | A Branch-and-Price-and-Cut Algorithm for Heterogeneous Pickup and Delivery Problems with Configurable Vehicle Capacity[END_REF]. In their problem, vehicles were not allowed to change configurations all along the route. This assumption grant configurations to be treated as vehicles types with some extra dimensioning constraints. DARP-RC instead, by allowing reconfigurations, introduces the challenge of tracking each configuration for every visited node. Vehicles can have one or several configurations. Each configuration is characterized by the capacity for each user type. Consider for example vehicle in Fig. 1.

In the first configuration it can handle 7 seated people and 1 wheelchair; in the second one, 6 seated people and 2 wheelchairs; in the third one, 4 seated people and 3 wheelchairs. Note that there is no linear relationship between the capacity in the two types of users (one wheelchair cannot be simply converted into one or two seat spaces). Also that unused chairs are folded and not removed from the vehicle.

Example.

The following example illustrates how vehicles with reconfigurable capacity can reduce operating cost. Consider a vehicle with two configurations {c 1 , c 2 } as shown in Fig. 2. The first configuration consists of 2 seats and 1 wheelchair space. The second configuration has 4 seats only. Users a and c go to destination M 1 while b, d, e and f go to M 2. In order to satisfy all user demands, the reconfigurable vehicle can follow the route

D → a → b → c → M 1 using configuration c 1 and d → e → f → M 2 → D with configuration c 2 .
Performing the same route without reconfiguring capacity would imply making an extra detour d → M 2 → d using configuration c 1 only (see dotted line); therefore increasing transportation costs.

Problem Description

The DARP-RC is defined as a network composed of a set V of vertices containing the set O + of vehicle starting depots, the set O -of vehicle arrival depots, the set P of pickup locations and the set D of delivery locations. Without loss of generality, we address the case of morning routes where the set P corresponds to people home or any nearby address, and the set D corresponds to MSEs.

The set of users is partitioned into several categories u ∈ U, i.e. seat, wheelchairs. An user request r ∈ R is composed of a pickup at location p r ∈ P, a delivery at location d r ∈ D and the number q ru of persons of each type u ∈ U to be transported. Moreover, each request has a maximum ride time Tr that guarantees a given quality of service.

The fleet of vehicles is homogeneous. Each vehicle has the set C of possible configurations that provides different capacity vectors. Hence, the vehicle capacity is expressed by quantities Q c u representing the maximal number of users type u ∈ U that can be transported at a time by the same vehicle when using configuration c ∈ C. Without loss of generality, we consider that if i ∈ P represents a pickup node, then the node i+n is the delivery node associated with the same request. It is assumed that different nodes can share the same geographical location.

Consider the directed graph G = (A, V) shown in Fig. 3. V = P ∪ D ∪ O corresponds to the set of nodes A is the set of arcs connecting the nodes. Every arc (i, j) in the graph represents the shortest path in time between nodes i and j. Its travel time is denoted as t ij while its length is denoted as d ij .

Note that reconfiguration time is not considered because it is negligible compared with service times and it can be performed at MSEs while massive drop-offs take place.

The objective function is to minimize the total transportation cost. Three main costs are considered: a fixed cost associated with each vehicle (amortization cost), a time related cost proportional to the total route duration (driver cost) and a distance related cost (vehicle fuel and maintenance costs).

RESOLUTION METHOD

An exact method has been proposed in [START_REF] Qu | A Branch-and-Price-and-Cut Algorithm for Heterogeneous Pickup and Delivery Problems with Configurable Vehicle Capacity[END_REF] to solve the heterogeneous PDP with configurable vehicle capacity for up to 50 requests. For large scale applications usually found in real-life situations, metaheuristics provide a good alternative. We propose a matheuristic that combines the metaheuristic LNS with SCP. It will be denoted as LNS-SCP.

Matheuristic framework (LNS-SCP)

The matheuristic framework consists of an LNS with a nested SCP solved periodically. LNS was first proposed by [START_REF] Shaw | Using constraint programming and local search methods to solve vehicle routing problems[END_REF] and introduced under the name ruin and recreate in [START_REF] Schrimpf | Record breaking optimization results using the ruin and recreate principle[END_REF]. In LNS the current solution is improved following an iterative process of destroying (i.e. removing parts of the solution) and repairing the current solution. This process is repeated until a stopping criterion is reach. In our case the stopping criterion is a maximal number of iterations or a maximum computational time.

The potential of the approach was revealed by Ropke and Pisiuger [START_REF] Ropke | An Adaptive Large Neighborhood Search Heuristic for the Pickup and Delivery Problem with Time Windows[END_REF] who proposed an ALNS consisting in multiple search operators adaptively selected according to their past performance in the algorithm.

Algorithm

s ← s If (cost of s is better than cost of s ) s * ← s If (it modulo η = 0) { /*set covering component*/ s ← solve set covering problem(Ω) Update s * ← s if s is cheaper than s* Update s ← s if s is cheaper than s Perform pool management } it ← it + 1 } Return s * -End -
LNS manages 3 solutions at a time: the current solution s, the temporary solution s generated after destroying and repairing a copy of s, and the best solution found so far s * . LNS consists in 3 fundamental steps:

(1) Determine the number of requests to remove Φ. We first randomly select the percentage of requests to be removed in the interval [α, β].

(2) Destroy and repair the current solution with a randomly selected operator σ -∈ Σ -and σ + ∈ Σ + respectively. This step results in a new temporary solution.

(3) Accept or reject the new temporary solution s , according to the record-to-record criterion of [START_REF] Dueck | New optimization heuristics[END_REF]: if objective(s ) ≤ (1+χ) * objective(s * ), where χ is a small positive value, s is accepted as the new current solution.

The SCP component is performed every η iterations. The purpose of the SCP component is to correct the LNS bias of discarding good routes that are part of costly solutions. Every new route is a candidate to be stored into a pool Ω of routes. Implementation details are presented in Section 4.3.

LNS operators

A key aspect of the LNS are the set of destroy and repair operators. With the goal of keeping the framework as simple as possible, we try to reduce the number of used operators without sacrificing the solution quality. At the end, only 2 destroy and 2 repair operators were kept in the framework based on the performance obtained in benchmark instances in Section 5.

Destroy operators determine the set of requests to be removed from the solution according to a certain criterion. In the framework of [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF] 7 destroy operators are proposed. After testing our framework on literature instances we find out that random removal and historical node-pair removal were sufficient to obtain competitive results. For details about the implementation of these operators, please refer to [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF].

Repair operators rebuild a partially destroyed solution in order to restore a complete feasible solution. This operation consists in reinserting nodes one by one from the request bank into the solution according to a specific criterion. Every insertion must satisfy all problem constraints. In our case time windows, ride times and capacity constraints have to be respected. If the operator does not fully repair the solution due to feasibility requirements, the objective function is strongly penalized by multiplying it by a big constant value. The two most common repair operators for the DARP are the cheapest insertion and the k-regret heuristics, both employed in the LNS-SCP. We implemented the k-regret heuristics with values of k varying from 2 to 4.

Set Covering Problem (SCP)

In LNS, a solution is rejected solely based on its cost with respect to the best solution so far. A rejected solution may contain some good routes, which are also removed. This issue is addressed by storing the routes found by LNS and using them in a SCP to find new solutions. In the following lines, we present the mathematical model and key components for its implementation.

Let Ω be the set of routes in the pool collected through the LNS iterations, and V ω ∈ R + the cost of route ω ∈ Ω. To describe the itinerary followed by each route, we define the values R rω , which are set at 1 if request r ∈ R is served by route ω ∈ Ω, and 0 otherwise.

The set covering problem aims at determining the value of binary variables y ω , where y ω = 1 if route ω ∈ Ω is part of the solution and 0 otherwise. The set covering problem is defined by the following model.

min ω∈Ω V ω y ω , (1) s
.t. ω∈Ω R rω y ω ≥ 1 ∀r ∈ R, (2)
y w ∈ {0, 1} ∀ω ∈ Ω. (3)
In every iteration, after calling the repair operator, the current routes are memorized in the pool Ω. In order to reduce the number of variables for the SCP, only nondominated routes are saved according to Proposition 1. Proposition 1. (Route dominance). Route ω 1 dominates route ω 2 , if ω 1 visits the same set of nodes as ω 2 at a lower cost.

The SCP is solved with a MILP solver every η iterations.

The solver is initialized with the best known solution and solved given a time limit T limit . This implies the SCP is not always solved to optimality. If the obtained solution is better that the current solution s then the current solution is updated. Otherwise s remains unchanged.

Similarly, the best solution s * is updated if a better solution is found. As constraint (2) of the set covering model allows request duplicates in the output solution, all duplicates are removed to obtain a consistent solution.

The cheapest one is always conserved.

If the solver fails to find an optimal solution within the time limit T limit , the pool is cleared and filled again with the routes of the best known solution. This step is refereed as pool management in Alg.1.

RESULTS

In order to determine the added value of the set covering component, we implemented three LNS variants: a classic LNS (using operators: k-regret, random removal, historical node-pair removal, worse removal, time related removal and distance related removal), an Adaptive LNS using the same set of operators and the proposed LNS-SCP (using only best insertion, k-regret, random removal and historical node-pair removal). In all cases k=2,3,4 regret was used. The MILP solver for solving the SCP is IBM Ilog Cplex 12.6, running on a single thread. The SCP is solved every η = 1000 iterations with a time limit of T limit = 3 seconds. The acceptance criteria is set to χ = 5% and the percentages used in destroy operators are α = 10% and β = 40% as in [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF].

In order to test the framework, a set of 9 instances of different sizes is proposed. The first number in the instance name refers to the number of considered requests. The service time for valid people was set to 2 minutes for pickup and 1 minute for delivery. While the service time for wheelchairs was 2 minutes for pickup and 5 minutes for delivery. There is not time window in pickup locations.

There is a single depot for vehicles with infinite capacity. Vehicles are all of the type illustrated in Fig 1.

Two experiments were completed. One limiting the number of iterations (Table 1) and the other one limiting the computation time (Table 2). In both cases it is compared LNS, ALNS, and LNS-SCP metaheuristics. 5 runs per instance type were considered. The average gap (Gap avg ) was computed with respect to the best known solution (BKS) found throughout both experiments.

In Table 1 can be observed that in most of cases, SCP component improves the objective function on average in 2.96% (= 3.99% -1.03%) with regard to the LNS and in 1.21% (= 2.24% -1.03%) with respect to the ALNS version of the algorithm. By limiting the computational May 2017

time to 30 minutes it can be observed a similar behavior as shown Table 2, however with a higher gap for the lasts instances. This is expected as the computation time is not big enough as to obtain good results in big instances.

Looking at very large instances, in both benchmarks, ALNS outperforms the other methods which may indicate that ALNS-SCP may be relevant for evaluation. On going experimentation also aim to establish the best value of parameter T limit . In particular a greater value for this parameter may be beneficial for large size instances. Table 3 shows detailed information of the best solution found for each instance. It shows computational time (t), the best cost among the five runs (obj), the number of reconfigurations (rec), the number of routes in the solution(routes) and the average number of requests per route (requests/route). The maximum number of reconfigurations was 1 for most of instances. These results is influenced by the vehicle fixed cost and maximum ride time constrains. Nevertheless, a deeper study should be done to determine precisely the key factor for reconfiguration. Finally we can observe the average number of requests per route of 7.31 that is higher than the vehicle capacity.

CONCLUSION

Through the paper we have described the dial-a-ride problem using vehicles with en-route reconfigurable capacity and a solution procedure based on LNS and SCP. By analyzing the performance of the LNS-SCP on real life instances, we could observe some significant gains when compared with LNS and ALNS.

We also observed that best LNS-SCP solutions reconfigure en-route at most once for most of the instances. The perspectives are therefore to extend the study with heterogeneous vehicles and to make an exhaustive evaluation to characterize the key factors of vehicle reconfiguration.

From the solution method side, we aim to establish the best value for the SCP time limit in particular for large size instances and include the ALNS-SCP for evaluation.

Fig. 1 .

 1 Fig. 1. Example of reconfigurable vehicle (source: www.handynamic.fr)

  Fig. 2. Comparison of routes with and without capacity reconfiguration 2

  1 shows the general structure of the LNS-SCP. Algorithm 1. (The LNS-SCP framework).Input: Σ -: set of destroy operators, Σ + : set of repair operators, η: nb. of iterations between two calls to the SCP. Output: best solution found s * .

	-Begin-
	Pool of routes: Ω ← ∅
	Request bank: B ← ∅
	s, s : current and temporary solutions
	it ← 0
	While (termination criterion not met) {
	s ← s
	Destroy quantity: randomly select a value Φ
	Operator selection: select σ

-∈ Σ -and σ + ∈ Σ + Destroy: apply σ -to remove Φ requests from s Copy the Φ requests into B Repair: apply σ + to insert requests from B into s Ω ← Ω ∪ routes(s ) /* add routes into the pool */ If (acceptance criterion is met)

Table 2 .

 2 Performance comparison of metaheuristics in 30 minutes over 5 runs

Table 3 .

 3 Best results