B.-H Nguyễn

J P F Trovão

The combination of batteries and supercapacitors is promising in electric vehicles context to minimize battery aging. Such a system needs an energy management strategy (EMS) that distributes energy in real-time for real driving cycles. Pontryagin's minimum principle (PMP) is widely used in adaptive forms to develop real-time optimization-based EMSs thanks to its analytical approach. This methodology leads to an off-line optimal solution which requires an extra adaptive mechanism for real-time applications. In this paper, a simplification of the PMP method is proposed to avoid the adaptation mechanism in real-time. This new EMS is compared to well-known conventional strategies by simulation. Furthermore, experimental results are provided to assess the real-time operation of the proposed EMS. Simulation and experimental results prove the advantages of the proposed approach by a reduction up to 50% of the batteries rms current on a real-world driving cycle compared to a battery-only EV.

Index Terms-Energetic macroscopic representation, Energy storage system, Hardware-in-the-loop, Optimal control. (EVs) are the future of our world due to limited fossil fuel sources and environmental issues. Energy storage, however, has a direct impact on the autonomy and the price of an EV. Batteries cost, energy density, and especially life-time are crucial for EV market. Adding supercapacitors (SCs), which forms hybrid energy storage systems (H-ESSs), is a promising solution [START_REF] Tie | A review of energy sources and energy management system in electric vehicles[END_REF]. Being known as a kind of high power-density sources, SCs can support batteries to reduce the degradation caused by current peaks.

NOMENCLATURE

Variables

The H-ESSs requires an energy management strategy (EMS) to distribute the energy in a fair way between batteries and SCs. EMS is an interesting topic attracting numerous efforts reported in the literature. They can be classified into rule-based and optimization-based strategies [START_REF] Salmasi | Control strategies for hybrid electric vehicles: Evolution, classification, comparison, and future trends[END_REF], [START_REF] Wirasingha | Classification and review of control strategies for plug-in hybrid electric vehicles[END_REF].

Rule-based EMSs can be fuzzy [START_REF] Zandi | Energy management of a fuel cell/supercapacitor/battery power source for electric vehicular applications[END_REF], [START_REF] Talla | Position-Based T-S Fuzzy Power Management for Tram With Energy Storage System[END_REF] or deterministic [START_REF] Wang | Optimal design and real-time control for energy management in electric vehicles[END_REF]- [START_REF] Machado | Effectiveness of Supercapacitors in Pure Electric Vehicle Using a Hybrid Meta-Heuristic Approach[END_REF], in which filtering strategy is among the most classical ones [START_REF] Thounthong | The benefits of hybridization[END_REF]- [START_REF] Song | Energy management strategies comparison for electric vehicles with hybrid energy storage system[END_REF]. Rule-based EMSs are often simple and easy to be implemented in real-time. Nevertheless, these methods require knowledge about the system behaviors and/or developers' expertise. That leads to non-optimal solutions. Real

E

On the other hand, optimization-based methods are more and more attractive due to their ability to minimize given cost functions. They are sorted as global optimal strategies and real-time sub-optimal ones [START_REF] Wirasingha | Classification and review of control strategies for plug-in hybrid electric vehicles[END_REF]. Global optimal EMSs are often developed by dynamic programming (DP) [START_REF] Koot | Energy Management Strategies for Vehicular Electric Power Systems[END_REF], Pontryagin's minimum principle (PMP) [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF]- [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF], or other optimization methods e.g. [START_REF] Hu | Longevityconscious dimensioning and power management of the hybrid energy storage system in a fuel cell hybrid electric bus[END_REF]. These strategies are optimal, however, only for off-line simulation due to the need of knowing the driving cycle in advance. Methods like DP require high computational effort which is hard for real-time application.

By contrast, real-time EMSs are sub-optimal but usable for real-time operation. Many optimization methods have been applied for developing EMSs such as model predictive control (MPC) [START_REF] Gomozov | Adaptive Energy Management System Based on a Real-Time Model Predictive Control with Non-Uniform Sampling Time for Multiple Energy Storage Electric Vehicle[END_REF], meta-heuristic methods [START_REF] Trovão | A multi-level energy management system for multi-source electric vehicles -An integrated rule-based meta-heuristic approach[END_REF], or adaptive PMP [START_REF] Kessels | Online Energy Management for Hybrid Electric Vehicles[END_REF]- [START_REF] Nguyen | Optimal control based algorithms for energy management of automotive power systems with battery/supercapacitor storage devices[END_REF]. They normally offer high performances for the managed systems. Their drawbacks are often on the complexity of the strategies. For instance, MPC and metaheuristic approaches need to do a lot of calculations to solve the optimization problems in every sampling time.

Adaptive PMP strategies, which are mostly -control [START_REF] Nguyen | Optimal control based algorithms for energy management of automotive power systems with battery/supercapacitor storage devices[END_REF], are less costly in terms of computational effort thanks to their analytical approach. The common way is to compute an openloop off-line optimal solution. A feedback of the state variable then leads to closed-loop adaptation of real-time evolutions. However, these methods require an important theoretical development based on a known-in-advance driving cycle that leads to sub-optimal solutions for real driving cycles not a priori known. A new approach has been proposed without adaptation mechanism [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF]. It consists in using a part of the PMP form including the state variable that indirectly induces a feedback (i.e. an adaptation). However, this method has been only studied in simulation; its real-time operation has not been demonstrated. Moreover, only a New European Driving Cycle (NEDC) driving cycle has been considered; its adaptation to other driving cycles has not been studied. Finally, no comparison with other real-time EMS was achieved.

The papers objective deals with an improvement and the positioning of the PMP-based EMS without adaptation mechanism proposed in [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF]. First, the co-state variable is determined in a new way. Second, the proposed EMS is compared with other classical EMS with different driving cycles to define its performances. Finally, its real-time operation is demonstrated by experiments using a laboratory set-up.

Section II describes the studied H-ESS, including modeling, control organization, and conventional EMSs. Section III presents the new EMS. Comparisons with classical EMSs are provided by simulation in Section IV. Experimental validations of the proposed EMS are presented in Section V.

II. CONVENTIONAL STRATEGIES OF THE STUDIED SYSTEM

A. Studied System Description 1) Configuration

The studied system is a semi-active battery/SCs H-ESS for EVs (Fig. 1). The batteries are coupled in parallel with a SCs subsystem. The SCs are connected in series with an inductor and a chopper. The H-ESS supplies the traction subsystem. Since this work focuses on the H-ESS, the traction subsystem is considered as an equivalent current source.

2) Modeling

All the mathematical equations for modeling are provided in Table I. The model representation can be found in the previous work [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF].

The batteries can be modeled by [START_REF] Tie | A review of energy sources and energy management system in electric vehicles[END_REF]. It outputs the batteries voltage ubat and inputs the batteries current ibat. There is a voltage drop in the equivalent series resistance (ESR) rbat. The ESR is in fact different in charging and discharging modes; and it is a non-linear function of the state-of-charge (SoC). The open-circuit voltage (OCV) ubat OC is also a non-linear function of the SoC SoCbat. In this work, the relationships rbat(SoCbat) and ubat OC(SoCbat) are given by look-up tables. The SoC is calculated by Coulomb counting with the capacitance Cbat.

Unlike battery, the SoC of SC is proportional to its voltage. Thus, the SCs model can be simply given by [START_REF] Salmasi | Control strategies for hybrid electric vehicles: Evolution, classification, comparison, and future trends[END_REF] with their voltage uSC, current iSC. The current iSC charges/discharges the SCs with the capacitance CSC. The SC resistance rSC causes an internal voltage drop and is quite independent of the SCs SoC. In this study, the rSC is given from the manufacture's datasheet.

The inductor, with its inductance L and internal resistance rL, is modeled by [START_REF] Wirasingha | Classification and review of control strategies for plug-in hybrid electric vehicles[END_REF]; in which uch sc is the chopper voltage.

The chopper model is expressed by [START_REF] Zandi | Energy management of a fuel cell/supercapacitor/battery power source for electric vehicular applications[END_REF]. The voltages ubat and uch sc are linked by the modulation function mch sc, as well as the current iSC and ich sc.

Equation ( 5) models the parallel connection, in which itrac is the traction current. The machine drive subsystem imposes itrac given by [START_REF] Wang | Optimal design and real-time control for energy management in electric vehicles[END_REF] which is the static model of the electrical machine drive. The efficiency ηm of the electrical machine drive is given by an efficiency map [START_REF] Dépature | Characterisation of the electric drive of EV: on-road versus offroad method[END_REF].

The mechanical transmission, combining final drive ratio and wheel transmission, is modeled by [START_REF] Carter | Optimizing for efficiency or battery life in a battery/supercapacitor electric vehicle[END_REF]. The vehicle chassis dynamics is addressed by the second Newton law of the relationship between the traction force, the resistant force, and the vehicle velocity given by [START_REF] Armenta | An advanced energy management system for controlling the ultracapacitor discharge and improving the electric vehicle range[END_REF]. The total mass Mveh is considered as the summation of the vehicle net weight and two passengers.

The environment imposes the total resistant force to the vehicle as a combination of rolling resistant force, air drag resistant force, and gravitational resistant force cause by slope as given in [START_REF] Machado | Effectiveness of Supercapacitors in Pure Electric Vehicle Using a Hybrid Meta-Heuristic Approach[END_REF].

3) Control Organization

The control scheme can be systematically conducted by the inversion of the system [START_REF] Bouscayrol | Graphic Formalisms for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF]. It is realized by equations given in the bottom part of Table I.

k TT T k T i T u                 (6) trac m trans m veh trans F T k vk      (7)   veh trac res 0 veh 1 t v F F dt M    (8)   2 res roll veh veh wind veh 0.5 x F k M g c A v v M g      (9)
Control ch sc ref trac meas bat ref i i i  (10) ch sc ch sc ref ch sc SC ref ch sc k i i m   (11)     ch sc ref SC meas P sc SC ref SC meas I sc SC ref SC meas 0 i t i u u k i i k i i dt      (12) ch sc ref ch sc bat meas u m u  (13)   trac ref res meas I veh ref veh meas P veh meas 0 t vv F F k v v dt k v      (14) trac ref m ref trans F T k  (15)
The coupling relationship can be directly inverted [START_REF] Thounthong | The benefits of hybridization[END_REF]. The chopper current reference ich sc ref is conducted from the batteries current reference ibat ref by compensating the measured traction current itrac meas.

The ich sc ref then enters to the direct inversion of the chopper current relationship which outputs the SCs current reference iSC ref as in [START_REF] Schaltz | Influence of battery/ultracapacitor energy-storage sizing on battery lifetime in a fuel cell hybrid electric vehicle[END_REF].

The indirect inversion of the inductor is achieved by a conventional proportional-integral (PI) current controller [START_REF] Allègre | Flexible real-time control of a hybrid energy storage system for electric vehicles[END_REF]; in which uch sc ref is the chopper voltage reference, uSC meas the SCs voltage measurement, iSC meas the SCs current measurement, kP i sc and kI i sc the PI coefficients, respectively. Finally, the modulation function mch sc is computed by [START_REF] Song | Energy management strategies comparison for electric vehicles with hybrid energy storage system[END_REF] with the measured batteries voltage ubat meas.

For control of vehicle dynamics, ( 8) is firstly inverted to obtain the velocity control loop. Since the chassis is an accumulation element which cannot be directly inverted, indirect inversion is required. It is realized by a closed-loop control [START_REF] Koot | Energy Management Strategies for Vehicular Electric Power Systems[END_REF] with the measurement of the vehicle velocity as the feedback. The electrical machine torque reference is then obtained by directly inverting the relationship [START_REF] Carter | Optimizing for efficiency or battery life in a battery/supercapacitor electric vehicle[END_REF] to obtain the reference torque [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF].

B. System Reduction for Energy Management

Energy management is at the higher layer than local control. Besides, different layers have different dynamics that need to be treated by different sampling time scales [START_REF] Trovão | A multi-level energy management system for multi-source electric vehicles -An integrated rule-based meta-heuristic approach[END_REF]. Furthermore, once an inner-loop subsystem is properly controlled, it can be considered to respond immediately to the reference. Hence, it is appropriate to reduce the fast dynamical subsystems by equivalent static model [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF], [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF].

The inductor, the chopper, and their control can be reduced to an equivalent conversion element. This equivalent element is considered without delay by assuming perfect response of the current control loop:

conv conv SC SC bat ch sc ch sc ch sc ref k u i u i ii         ; ( 16
)
in which ηconv is the efficiency of the DC/DC converter, kconv the coefficient depending on the power direction. This reduced is organized by using Energetic Macroscopic Representation (EMR) (Fig. 2). EMR is a graphical formalism for control organization of energetic systems [START_REF] Bouscayrol | Graphic Formalisms for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF]. It has been widely used for control of electric and hybrid vehicles [START_REF] Lhomme | Energy Savings of a Hybrid Truck using a Ravigneaux Gear Train[END_REF].

C. Filtering Strategy as a Traditional Rule-Based EMS

Filtering strategy is developed based on the fact of different frequency characteristics of energy storages [START_REF] Thounthong | The benefits of hybridization[END_REF]. It can be realized by using a low-pass filter (LPF) as follows:

bat ref trac meas LPF 1 1 ii s    ; ( 17
)
where τLPF is the time constant of the LPF given by

LPF LPF 1 2 f    ; ( 18
)
in which fLPF is the cut-off frequency of the LPF. The frequency can be determined by considering Ragone plot [START_REF] Allègre | Flexible real-time control of a hybrid energy storage system for electric vehicles[END_REF], by adaptive studies [START_REF] Florescu | Adaptive frequency-separation-based energy management system for electric vehicles[END_REF], or by iterative examines [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF].

D. -control as an Adaptive PMP-based Strategy

In this paper, for comparison purpose, the -control developed in [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF] is addressed due to the same studied system configuration. In that work, the SCs energy is considered as the state variable for the reduced model. All the variables are then transformed into the batteries side. As a result, the obtained control law is conducted as follows:

    ch sc bat OC bat trac meas bat ref bat bat OC ch sc 0 21 k k u r i i ru                 ; ( 19
)
where  is the co-state variable, 0 the initial co-state which is pre-computed by iterative simulations, Δ the adaptation. The adaptation is conducted by a classical PI controller to control the squared SCs voltage uSC 2 .

E. Dynamic Programming as an Off-line Optimal Benchmark

Thanks to the advantage of obtaining a global optimal solution, DP is used to deduce a benchmark for comparison. The method is based on solving the Bellman equation:

              * , * 1, , min , D kN uk kN g x k w k J x k J f x k w k               ; ( 20
)
in which J is the cost function to be minimized, N the total number of time steps, f the dynamical function of the system, gD the current accumulative cost, x(k) and w(k) the state and control variables at the k th step, respectively.

In the right-hand side, the term gD is the cost-to-go from the current k th step to the next (k+1) th step. The term Jk+1, N * is the known-in-advance optimal cost-to-go from the (k+1) th step to the final step N th . DP is therefore a backward method that can be computed off-line only.

III. PROPOSED REAL-TIME STRATEGY

A. Problem Formulation 1) Mathematical System Dynamics

To apply optimal control theory, the energy management problem should be formulated as an optimal control problem. For that, the system dynamical model serves as a set of equality constraints. The limitations of the state and control variables are inequality constraints.

In this work, the limitations are treated by the conventional method addressed in [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF]. Thus, only the equality constraint, i.e. the dynamical model of the system, needs to be defined.

From the reduced representation (see Fig. 2), the reduced mathematical model of the system can be deduced as follows:

conv ch sc ref SC bat SC SC k i d uu dt C u   ; ( 21
)
where uSC is the state variable, ich sc ref the control variable, ubat the disturbance input. Here, the resistance effect is neglected.

2) Performance Criterion -Cost Function

Since the objective is to reduce the batteries degradation, the rms value of its current is a relevant performance index:

f 0 2 bat rms bat f0 1 t t I i dt tt    ; ( 22
)
in which t0 and tf are the initial and the final time, respectively. Furthermore, quadratic functions are of interest for many optimization techniques [START_REF] Bryson | Applied Optimal Control: Optimization, Estimation, and Control[END_REF]. Hence, the cost function of the studied energy management problem can be defined by:

f o 2 bat t t J i dt   . ( 23
)

B. Strategy Development Using PMP 1) Methodology

First, the general background of PMP is briefly described. For a dynamical system  

,,

d x f x w t dt  ; ( 24
)
where x is the state variable and w control variable; with the cost function given by

    f o ff , , , t t J h x t t g x w t dt     ; ( 25
)
in which h reflects the cost of final state x(tf) and g the integral cost of the progress; a Hamiltonian function can be defined by

    , , , , H g x w t f x w t   ; ( 26
)
where  is the co-state variable.

The PMP states that if a control law w * is optimal, it must satisfy the following three necessary conditions:

  * * * * , , , H x w t d x dt      ; ( 27
)   * * * * , , , H x w t d dt x      ; ( 28
)     * * * * * , , , , , , H x w t H x w t   . ( 29
)
Conditions ( 27) and ( 28) compose the dynamics of a Hamiltonian system. Condition [START_REF] Florescu | Adaptive frequency-separation-based energy management system for electric vehicles[END_REF] means the minimization of the Hamiltonian function. Conventional methods satisfy all conditions to obtain an off-line open-loop control law; an additional state feedback is then used for adaptation in realtime.

An alternative approach is proposed in [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF]. Only condition (29) is used to directly deduce a real-time strategy. The Hamiltonian function contains a relationship between the control variable ich sc ref and the state variable uSC. Hence, applying (29) by a partial derivative with respect to the control variable can eventually carry out an analytical closed-form strategy. The co-state variable  is then determined based on its physical meaning without considering the Hamiltonian system dynamics ( 27) and [START_REF] Lhomme | Energy Savings of a Hybrid Truck using a Ravigneaux Gear Train[END_REF]. The obtained strategy is indeed sub-optimal as the PMP theory is not strictly satisfied.

2) Strategy Development

From [START_REF] Trovão | A multi-level energy management system for multi-source electric vehicles -An integrated rule-based meta-heuristic approach[END_REF], assuming an ideal converter, the Hamiltonian is defined by:

2 bat ch sc ref bat SC SC ui Hi Cu   . ( 30
)
Condition ( 29) is realized by

ch sc ref 0 H i    . ( 31
)
The calculation of the left-hand side of ( 31) is given from: bat trac ch sc i i i  .

(32) Besides, itrac is the disturbance without dependence on the control variable ich sc ref.

With assumption [START_REF] Li | Design and analysis of a fuel cell supercapacitor hybrid construction vehicle[END_REF], that leads to

  2 bat bat ch sc ref 2 i i i    . ( 33
)
The partial derivative to be calculated therefore becomes

SC ch sc ref SC bat ch sc ref bat 2 SC ch sc ref SC 2 u iu ui H i C iu          . ( 34
)
The partial derivative of uSC with respect to ich sc ref is calculated by:

SC ch sc ref SC SC SC SC SC bat bat SC bat SC 11 . u i i u i u i u u u u u           (35)
The partial derivative of iSC with respect to uSC is adopted from [START_REF] Nguyen | Optimal control based algorithms for energy management of automotive power systems with battery/supercapacitor storage devices[END_REF] as follows:

SC SC 2 SC SC SC SC 4 ii u u r P     . ( 36
)
Applying (36) to [START_REF] Allègre | Reduced-scale-power hardware-in-the-loop simulation of an innovative subway[END_REF], that leads to: (39) However, considering the control scheme of Fig. 2, the strategy delivers the batteries current reference ibat ref. Then, ich sc ref is obtained by compensation of the measured traction current as in [START_REF] Thounthong | The benefits of hybridization[END_REF]. Thus, the control law should be:

u u r u i u u r u i u i u i C u      . ( 42
)
Once the co-state variable  is determined, (42) will be ready to serve as a real-time strategy.

3) Co-state Variable  Physical Meaning and Determination

This approach determines  not by the Hamiltonian dynamics ( 27) and ( 28), but by its physical meaning. The physical essence of the co-state variable has been pointed out in the literature [START_REF] Ross | A Primer on Pontryagin's Principle in Optimal Control[END_REF]. It is an equivalent factor to convert the dynamical function of the system to an equivalent cost to be minimized. Such a physical explanation has been also figured out by previous works, e.g. [START_REF] Nguyen | Optimal control based algorithms for energy management of automotive power systems with battery/supercapacitor storage devices[END_REF]. It is used here to calculate .

From the Hamiltonian [START_REF] Bryson | Applied Optimal Control: Optimization, Estimation, and Control[END_REF], whereas the cost term ibat 2 has a dimension of Ampere squared, the dimension of the dynamics term is Ampere per Farad. To perform a proper subtraction, the co-state variable, as an equivalent factor, should have a form of  = Capacitance  Current. On the right-hand side, the former term can be chosen as the SCs capacitance which is the only capacitive quantity of the studied system.

The latter term of current is more flexible to be determined. In the previous work [START_REF] Nguyễn | An Optimal Control-Based Strategy for Energy Management of Electric Vehicles using Battery/Supercapacitor[END_REF], it is given by the expected maximal batteries current when the SCs can no longer support the batteries. This value is set by the proportions of the batteries C-rate regarding the driving modes. This approach leads to a high-performance strategy validated by the reported simulation results. However, C-rate is an indicator of batteries energy rather than their power. The justification of this approach is therefore not so convinced.

To overcome the issue mentioned above,  is determined from the traction rms current:

SC trac rms Ci   ; ( 43
) in which, itrac rms depends on the driving modes and should be pre-calculated. The pre-calculation can be done with a standard driving cycle. In this work, the Worldwide harmonized Light duty driving Test Cycle (WLTC) is used. Based on the studied vehicle, the WLTC class 2 (Fig. 3) is used for the pre-calculation of the itrac rms.

Low

Middle High Fig. 3. WLTC class 2 used for co-state variable  determination.

4) Strategy Implementation

The strategy is in an analytical form with simple mathematical calculations. It is therefore easy to be implemented in real-time (Fig. 4). Three sensors are required for the measurements of the SCs voltage uSC, batteries voltage ubat, and SCs current iSC to compute the control variable ibat ref cal from (42) and (43). A supplementary sensor for itrac is used only when the SCs voltage reaches limitations. In this case, the ibat ref is switched from ibat ref cal to the itrac meas.

IV. COMPARATIVE EVALUATION BY SIMULATION

A. Simulation Conditions

The reference vehicle is Tazzari Zero (main parameters in Table II). This EV is originally supplied by a battery pack of Lithium-iron phosphate LiFePO4 (LFP) cells. The battery has been sized for a real driving range around 100 km. For studies of H-ESSs, it is modified to the semi-active configuration by adding a SC subsystem. SC and battery size are constant during the study to fairly compare the different EMS strategies.

Three standard driving cycles are considered for testing:  WLTC class 2, which is utilized to determine the co-state variable , with urban, rural, and highway parts. This is the result of a statistical study that addressed various standard and real-world driving cycles from countries in Europe, Asia, and America [START_REF] Tutuianu | Development of a World-wide Worldwide harmonized Light duty driving Test Cycle[END_REF].  NEDC with urban and highway parts. This driving cycle is commonly used in many studies because it can examine clearly the acceleration and deceleration behaviors of vehicles. This cycle is however too smooth in comparison with the real-world-based cycles.  ARTEMIS urban which is an urban cycle. This cycle addresses better than NEDC the fluctuations especially in urban areas. Moreover, a real-world driving cycle was recorded through driving around the campus of University of Lille. It can be considered as a rural cycle.

The evaluation criterion is the batteries rms current; the lower the ibat rms is, the better the battery life-time is. It is because the rms current cause self-heating in batteries due to Joule losses in their internal resistance. Whereas the degradation of the batteries increases with the rise of temperature [START_REF] Baghdadi | Lithium battery aging model based on Dakin's degradation approach[END_REF], [START_REF] Waldmann | Temperature dependent ageing mechanisms in Lithium-ion batteries -A Post-Mortem study[END_REF].

Five cases are compared: (i) conventional EV with battery only, (ii) filtering strategy as a popular rule-based EMS, (iii) -control as an adaptive PMPbased EMS, (iv) the proposed sub-optimal strategy, and (v) DP as an optimal benchmark.

Considering a previous work in the field [START_REF] Castaings | Comparison of energy management strategies of a battery/supercapacitors system for electric vehicle under real-time constraints[END_REF], this comparison uses the cut-off frequency fLPF = 48 mHz and the same -control parameters for adaptive PMP-based EMS. The simulations are carried out in MATLAB ® /Simulink ® .

B. Results and Discussions

Simulation results for comparative evaluation are reported in Fig. 5. The pure battery EV always leads to the highest battery rms current; meanwhile, the DP offers the lowest ibat rms thanks to its optimal solution. The H-ESS has significant effectiveness with the highly fluctuated driving conditions such as ARTEMIS urban and especially the studied real-world cycle. For the last one, the proposed strategy can reduce up to 50% of the batteries rms current that means life-time extension.

In the less fluctuating cycles like NEDC and WLTC, the H-ESS has a lower impact. In the case of NEDC, the ideal solution from DP can only save 11.4% of the ibat rms; whereas it is 10% by using the proposed real-time strategy. H-ESS seems more suitable for city cars than for the other sorts of vehicles working with smoother conditions.

The most important is that the proposed EMS leads to the best results for any tested cycle.

V. EXPERIMENTAL VALIDATION OF THE PROPOSED STRATEGY

A. Experimental Setup

Experimental validations are carried out by using reducedscale power hardware-in-the-loop (HIL) simulation. This study focuses on the energy management strategy of the battery/SC H-ESS. The traction subsystem therefore plays the role of generating the traction current which is the disturbance of the studied system. Thus, the traction subsystem, including the vehicle dynamics, can be emulated by using a controllable current source. Similarly, a controllable voltage source is used to emulate the batteries.

Fig. 6 illustrates the principle of the reduced-scale power HIL simulation used in this work. The traction model, including the vehicle inertia, generates the traction current reference. Via the power adaptation block, this reference is imposed to the reduced-scale experimental test bench. If the current of the emulator follows the generated reference, the emulator can emulate well the traction dynamics. It is similar for the batteries emulation.

Voltage and currents controllers are realized by using the classical PI controllers. The models and the controllers are implemented in the controller board (dSPACE DS1103). The current and voltage sources for emulations are realized by using SCs, inductors, a DC bus capacitor, and power electronics bidirectional choppers (see Fig. 7). The SCs are used to emulate the batteries and the traction subsystems thanks to their fast dynamics. It is noteworthy that the dynamics of the components used for emulation should be faster than that of the emulated system. The experimental setup is shown in Fig. 8. Main parameters of the components are given in Table III.

The voltages and the currents of all the SC banks and the emulated batteries voltage (on the capacitor Cbat emu) are measured by sensors. The currents ibat HIL, itrac HIL, and ich sc HIL cannot be directly measured because they are pulses as a result of choppers switching. They are estimated by using the corresponding measured SC currents iSC bat, iSC trac, and iSC. The estimation is obtained by using the average conversion relationship of the choppers (similar to (4)).

The power and the energy of the studied system are scaled down by a reduction ratio for power HIL implementation [START_REF] Allègre | Reduced-scale-power hardware-in-the-loop simulation of an innovative subway[END_REF]. The adaptation ratio is determined regarding the power and energy limitation of the available power components. It must ensure that the SCs used for traction emulation can absorb all the energy delivered by the SCs for the batteries emulation. The required energy storage capabilities are calculated with the WLTC, which is the longest driving cycle under study. Regarding that, the adaptation ratio is set to 35.

The emulated batteries voltage is equal to that voltage of the full-scale system. Thus, the currents itrac, ibat, and ich sc are scaled down exactly by the power ratio.

It is slightly more complicated for scaling the SC voltage. The SC energy must be scaled down by the same reduction ratio as that of the power. The simplest way is to reduce the SC capacitance by that ratio. The SC voltage range can be therefore remained as in the full-scale system. Unfortunately, the SC pack available in the laboratory has higher capacitance than it should do. Thus, the SC voltage range must be reduced. The upper boundary uSC max are set the same for both full-scale simulations and reduced-scale experiments. But the lower boundary uSC min of the reduced-scale system is higher than that of the full-scale system.

B. Results and Discussions

To validate the real-time performances of the proposed strategy, the real-world driving cycle is under study (Fig. 9).

The vehicle velocity of the studied driving cycle is imposed (Fig. 9). The simulation results of the full-scale studied vehicle present the evolutions of the batteries voltage, the SCs voltage, and the H-ESS currents (Fig. 10). The experimental results of the reduced-scale power HIL system present the corresponding evolutions (Fig. 11).

Regarding the power and energy reductions, it can be pointed out that the experimental results match simulation results. There is only a very small difference in the SCs voltage uSC. It is due to the difference between the internal resistance of the reduced-scale SC (see Table III) and that of the full-scale SC (see Table II). The fine-matched results verify the ability of the reduced-scale power HIL simulation in examining the real-time performances of the EMS.

The driving condition of the studied real-world cycle fluctuates very much; that is suitable for SCs. That is why they can perfectly protect the batteries from peak currents required by the traction subsystem. A reduction of about 70% of the batteries peak currents is reported in comparison to the demanded traction current.

The accuracy of the emulations can be estimated by

where x generally denotes itrac and ubat. The relative errors (see Table IV) are very small that confirms the controlled current and voltage sources can emulate well the traction subsystem and the batteries, respectively. Table IV also gives quantitative results to verify the performance of the proposed EMS in comparison between the full-scale off-line simulation and the reduced-scale HIL simulation. The batteries rms currents are reduced in comparison with the traction ones thanks to the H-ESS managed by the proposed strategy. The reductions are 49.6% and 49.8% for full-scale and reduced-scale systems, respectively. That verifies the reduced-scale HIL system examines exactly the behaviors of the full-scale studied system. 11. Experimental results of the reduced-scale power HIL system.

VI. CONCLUSION

This paper proposes and validates a new approach of using PMP to develop real-time EMS for H-ESSs in EVs. By using only, the third necessary condition of PMP, a closed-form solution containing the state variable can be deduced. Thus, no additional adaptation of the co-state variable is required for real-time applications.

This novel strategy has been compared to the conventional EMSs including filtering, -control, and DP. Simulation results have figured out that the proposed EMS gives the highest benefit in term of batteries rms current reduction. By comparing to the off-line optimal benchmark given by DP, the EMS has been verified to be close-to-optimal. Moreover, the real-time performances of the proposed strategy have been demonstrated by reduced-scale power HIL experiments.

Fig. 1 .

 1 Fig. 1. Studied system: a semi-active battery/SC H-ESS for EVs.

Fig. 2 .

 2 Fig. 2. Reduced EMR for energy management study.

Fig. 4 .

 4 Fig. 4. Implementation of the proposed EMS.

Fig. 5 .

 5 Fig. 5. Comparative evaluation of the studied EMSs.

Fig. 6 .

 6 Fig.6. Reduced-scale power HIL system for experimental validation.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Experimental system hardware configuration.

Fig. 9 .

 9 Fig. 9. Velocity of the studied real-world driving cycle.

Fig. 10 .

 10 Fig. 10. Simulation results of the full-scale vehicle.

  -Time Energy Management of Battery/Supercapacitor Electric Vehicles Based on an Adaptation of Pontryagin's Minimum Principle Bảo-Huy Nguyễn, Ronan German, João Pedro F. Trovão, Senior Member, IEEE, and Alain Bouscayrol, Member, IEEE

TABLE III REDUCED

 III -SCALE POWER HIL SYSTEM PARAMETERS

	Parameters		Values
	Adaptation ratio		35
	SCs for battery and traction emulators
	SC capacitance	CSC bat/trac	130 F
	SC nominal voltage	uSC bat/trac nom	54 V
	SC subsystem
	SC capacitance	CSC	14.5 F
	SC nominal voltage	uSC nom	60 V
	SC internal resistance	rSC	73.0 mΩ
	Inductor of the battery emulator
	Inductor inductance	Lbat emu	0.758 mH
	Inductor parasitic resistance	rL bat emu	0.20 Ω
	Inductor of the traction emulator
	Inductor inductance	Ltrac emu	0.752 mH
	Inductor parasitic resistance	rL trac emu	0.20 Ω
	Inductor
	Inductor inductance	L	0.751 mH
	Inductor parasitic resistance	rL	0.18 Ω
	Capacitor of the battery emulator
	Capacitor capacitance	Cbat emu	2200 uF

This work was supported by the Canada Research Chairs Program, the National Sciences and Engineering Research Council of Canada, and also with the support of the "Service de Coopération et d'Action Culturelle du Consulat Général de France à Québec" within the Samuel-De Champlain Program. This work has been achieved within the framework of CE2I project co-financed by European Union with the financial support of European Regional Development Fund (ERDF) and the French Region of Hauts-de-France. The review of this paper was coordinated by Dr. Mohamed Benbouzid.


Frank Pennekamp
email: frank.pennekamp@ieu.uzh.ch

Mikael Pontarp
email: mikael.pontarp@ieu.uzh.ch

Andrea Tabi
email: andrea.tabi@ieu.uzh.ch

Florian Altermatt
email: florian.altermatt@eawag.ch

Roman Alther
email: roman.alther@eawag.ch

Yves Choffat
email: yves.choffat@ieu.uzh.ch

Emanuel A Fronhofer
email: emanuel.fronhofer@umontpellier.fr

Pravin Ganesanandamoorthy
email: pravin.ganesanandamoorthy@eawag.ch

Aurélie Garnier
email: aurelie.garnier.research@gmail.com

Jason Griffiths

Suzanne Greene
email: e-segreene@mit.edu

Katherine Horgan
email: katherine.horgan@ieu.uzh.ch

Thomas M Massie
email: thomas.massie@outlook.com

Gian-Marco Palamara
email: gianmarco.palamara@gmail.com

Mathew Seymour
email: d-m.seymour@bangor.ac.uk

Biodiversity increases and decreases ecosystem stability

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

We performed a factorial manipulation of the diversity and composition of competing species (1 to 6 species, 53 unique community compositions) and temperature (six constant levels, modelled as a linear predictor) in microbial communities of bacterial consumers, and recorded community biomass dynamics over time. For each replicate we then calculated two stability components: resistance (= [total biomass at T˚C -total biomass at 15˚C] / [T˚C -15˚C] where T is the temperature of the replicate) and the temporal stability of biomass (inverse of coefficient of variation of community biomass). While these stability indices are widely used by empiricists, they should not be mistaken for mathematical definitions such as asymptotic resilience, which are more precise but also more restrictive [START_REF] Arnoldi | Resilience, reactivity and variability: A mathematical comparison of ecological stability measures[END_REF] .

Increased species richness caused greater temporal stability of total biomass (figure 1-a1) (linear mixed model of log inverse CV: richness effect size 0.33 with a standard error of 0.065) at all temperatures (extended data figure 1). Total biomass increased during the first week of the experiment and then declined over the next five weeks and total biomass was higher in more species-rich communities (figure 1-b, 1-c2, extended data table 1) (effect size for log richness 0.05 [units of mg/mL/log(species richness) unit] with 0.0096 standard error).

In contrast, increased species richness decreased resistance of total biomass to warming (figure 1-a2) (negative effect of log richness in a linear model, effect size of -0.006 [mg/°C/ log(species richness) unit] with a standard error of 0.0018). Richness negatively affected resistance measured on both absolute and relative scales (extended data figure 2). This effect was corroborated in analyses of total biomass by a negative interaction term between temperature and richness, which persisted through the experiment except during the first days (figure 1-c3) (log(richness) x temperature interaction of -0.0053 [units of mg/mL/°C/log(species richness) unit] with standard error of 0.00051) despite large variation in dynamics of total biomass (figure 1-b). This negative interaction reflects a stronger negative effect of temperature on total biomass (i.e. lower resistance) in richer communities (i.e. a richness-dependent response of total biomass to temperature).

Hence, temporal stability and resistance were negatively correlated across the species richness gradient (figure 1-d, RMA analysis with slope = -0.009, 95% CI = -0.0178 to -0.0051). Niche complementarity, statistical averaging, low overall response diversity, and possibly lower response diversity in more diverse communities were likely causes of the opposite effects of richness on temporal stability (extended data figure 3). The two stability components were, however, positively correlated within any single level of species richness (figure 1-d, extended data table 2). That is, composition variation without changes in species richness resulted in positively covarying temporal stability and resistance.

Next, we examined studies (including our own) measuring multiple stability components across diversity gradients based on a review by Donohue et al. (2016) [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF] (figure 2, extended data table 3 & 4). Seven of 30 comparisons show positive covariance, twenty show no covariance, and three showed negative covariance. Our study for the first time identifies negative covariation between resistance and temporal variability caused by intrinsic dynamics only. Although infrequently reported, negative covariation is disproportionately important because it complicates conclusions about and practical implications of effects of diversity on stability. Furthermore, these studies may be unrepresentative of the true prevalence of negative covariation, due to it being overlooked, publication bias towards positive diversity-stability relationships [START_REF] Ives | Stability and Diversity of Ecosystems[END_REF] or if the scale of analysis masks such covariation, e.g. within richness versus across richness.

A general mechanistic understanding of why different studies find different correlations would be a major step forward. Of the 30 pairs of stability components, only seven were accompanied by quantitative analyses of mechanism for both diversity-stability relationships (extended data table 4). Response diversity was implicated in five of these seven. Indeed, response diversity has been identified as an important driver of the resilience of ecological systems [START_REF] Suding | Scaling environmental change through the community-level: a traitbased response-and-effect framework for plants[END_REF]24 , and correlation among effect (i.e. high biomass production) and response (e.g. response to an environmental driver) traits, at least in the absence of strong interspecific interactions, might predict covariance among stability components [START_REF] Suding | Scaling environmental change through the community-level: a traitbased response-and-effect framework for plants[END_REF] . None of the studies involved manipulation of mechanism.

Negative covariation also raises the potential for non-monotonic effects of diversity on overall ecosystem stability. Ecosystem multifunctionality is the simultaneous consideration of several ecosystem functions. Similarly, overall ecosystem stability considers multiple stability components simultaneously [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF] (see extended data figure 4 for clarification of terms and concepts). We first calculated overall ecosystem stability as the sum of standardised individual ecosystem stability components [START_REF] Maestre | Plant species richness and ecosystem multifunctionality in global drylands[END_REF] . For the results of our experiment, in which one stability component increases linearly with diversity, and another decreases linearly (figure 1-a1,2), overall ecosystem stability is invariant with species richness since the two standardized components perfectly balance each other.

We further explored the implications of the negative covariation among stability components. Existing methods for calculating ecosystem multifunctionality include a mathematical function for converting values of each ecosystem function into a common currency. Following others [START_REF] Slade | The importance of species identity and interactions for multifunctionality depends on how ecosystem functions are valued[END_REF] , we assumed a logistic shape conversion function, and converted observed stability components into the common currency, with parameter values chosen to ensure that the range of the data includes the lower and upper asymptote. Here we present the influence of the value of the intercept (parameter Q) of the conversion function (figure 3a), which is similar to varying the threshold value when calculating ecosystem multifunctionality with the threshold approach [START_REF] Gamfeldt | Revisiting the biodiversity-ecosystem multifunctionality relationship[END_REF] .

The opposing stability-diversity relationships from our experiment, coupled with logistic conversion functions produced a hump-shaped relationship between overall ecosystem stability and diversity at low values of Q, a flat relationship at intermediate values, and a U-shaped relationship at higher values (figure 3b,c). Effects of richness on overall ecosystem stability were weaker at the replicate level than the richness level (figure 3b) in part due to relatively large variability among the replicates within richness levels and also due to positive covariance of temporal stability and resistance within diversity levels (figure 1d).

Finding hump-and U-shaped diversity-stability relationships fundamentally alters the importance of biodiversity, from it being unimportant (variation in biodiversity has no consequence), to being state dependent. The hump-shaped and U-shaped relationships occurred because measured stability components (e.g. resistance) cross the threshold of the conversion function at different levels of species richness (see link to interactive figure in Methods). Thus, negatively covarying stability components and non-linear conversion functions can produce a non-monotonic effect of diversity on overall ecosystem stability. Furthermore, we expect the same to apply to ecosystem multifunctionality when individual ecosystem functions negatively covary along a diversity gradient.

The chosen mathematical conversion function and its parameter values are therefore important. Ecosystem multifunctionality-diversity relationships can be negative, neutral or positive depending on the chosen threshold value [START_REF] Gamfeldt | Revisiting the biodiversity-ecosystem multifunctionality relationship[END_REF] . Similarly, different weightings of ecosystem functions can alter which species and interactions are considered important for multifunctionality [START_REF] Slade | The importance of species identity and interactions for multifunctionality depends on how ecosystem functions are valued[END_REF] . Other types and parameterisations of conversion functions might conceivably result in other diversity-stability relationships. For instance, weighting one ecosystem function (or stability component) highly, would lead to a multifunctionality / overall ecosystem stability measure that is essentially univariate. We view the choice and parameterization of conversion functions as an opportunity to tailor estimates of ecosystem multifunctionality and overall ecosystem stability to their policy and decision making applications [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF][START_REF] Slade | The importance of species identity and interactions for multifunctionality depends on how ecosystem functions are valued[END_REF][START_REF] Manning | Redefining Ecosystem Multifunctionality[END_REF][START_REF] Armsworth | The economic value of ecological stability[END_REF] . These conversion functions and estimates of overall ecosystem stability can then transform the perceived importance of biodiversity change for the sustainable delivery of multiple ecosystem services and may assist in translating the results of scientific studies into actionable information.

Methods section

Experimental methods

The experiment was conceived and designed to research biodiversity-ecosystem relationships, with a focus on questions on environmental gradients and effects on temporal changes / stability, including the one addressed in this article. We factorially manipulated temperature (15, 17, 19, 21, 23 and 25 °C) and species richness (1 to 6 species of bacterivorous ciliates: Colpidium striatum, Dexiostoma campylum, Loxocephalus sp., Paramecium caudatum, Spirostomum teres, and Tetrahymena thermophila) in 690 microcosms (250 mL Duran bottles). There were two incubators at each temperature. Manipulating temperature with a replicated gradient is recommended to harness the power of a regression design, while still allowing us to test for a nonlinear temperature effect [START_REF] Cottingham | Knowing when to draw the line: designing more informative ecological experiments[END_REF] . During prior testing the temperature of liquid in similar microcosms varied around the set-point temperature by 0.1˚C. Long-term ciliate cultures were kept at 15 °C, and warming usually decreases their carrying capacities but increases growth rates [START_REF] Leary | Testing a biological mechanism of the insurance hypothesis in experimental aquatic communities[END_REF] . For Colpidium, temperatures above 20 °C represent a stress that decreases population growth rate, whereas the other species are still within their tolerated thermal range.

To start the experiment, ciliates were grown to their respective carrying capacities at 20°C in bottles containing 1L of bacterized medium. Throughout the experiment, medium consisted of protist pellets (Carolina Biological Supplies, Burlington, NC, USA) at a concentration of 0.055g L -1 of Chalkley's medium in which the bacterium Serratia fonticola was grown as the resource for all ciliate consumers. Two autoclaved wheat seeds were added to each bottle for slow nutrient release [START_REF] Altermatt | Big answers from small worlds: a user's guide for protist microcosms as a model system in ecology and evolution[END_REF] .

Monocultures (species richness = 1) were initiated at a density of 3 individuals mL -1 in a total of 100mL medium. Multispecies communities (containing 2 to 6 species) were initiated with a total of 40 mL ciliate culture topped up with 60 mL fresh medium (100 mL culture in total). The 40mL cultures were assembled by adding a fixed fraction (i.e. 20 mL for 2 species, 13.33 mL for three species etc.) of each species at their specific carrying capacity, adopting a substitutive design. Since the number of possible species compositions exceeded the number of feasible experimental units, we used all possible compositions only for the monocultures, two and six species communities. For all other levels, species compositions were selected randomly from the set of all possible compositions such that all species occurred the same number of times, resulting in a total of 53 different compositions. Each level of species richness and composition was replicated at least twice, including an additional replicate for the two and three species level, and 3 additional replicates for the six species community resulting in 115 experimental units per temperature.

We sampled each experimental unit every day for the first 7 days, then 3 times per week for the following 50 days and a final sampling 7 days later, resulting in time series of 27 time points over a 57-days period. We used video sampling techniques to count and measure individual ciliates in all communities [START_REF] Pennekamp | Dynamic species classification of microorganisms across time, abiotic and biotic environments-A sliding window approach[END_REF] . For sampling, microcosms were taken out of the incubator, gently stirred to homogenize the culture and a sample was pipetted into a counting chamber. The counting chamber was covered with a lid and a 5s long video was taken under the microscope. The videos were subsequently processed with the R package BEMOVI [START_REF] Pennekamp | software for extracting behavior and morphology from videos, illustrated with analyses of microbes[END_REF] .

We derived community biomass by summing the biovolume of all individuals of a given species in a given community and multiplying biovolume with a constant density equal to water (i.e. 1g/cm 3 ). For each community, this resulted in a time series of community biomass. To avoid analysing monoculture time series whilst starting at a low fraction of the carrying capacity (3 individuals mL -1 ), compared to our multi-species culture species which started between 7 and 20 % carrying capacity (depending on richness), we aligned the monoculture time series such that they started at the day where they crossed 20% of their expected carrying capacity. Results regarding the effect of richness on temporal stability are qualitatively robust to exclusion of monoculture data (extended data figure 1c).

The temporal stability of each replicate community was calculated as the inverse of the coefficient of temporal variation in community biomass (mean community biomass / standard deviation of community biomass). Resistance was measured for each replicate as the absolute difference between total biomass in the replicate and the mean total biomass at 15˚C, divided by the temperature difference, with dimensions of milligrams (per ml) per degree Celsius. Resistance values close to 0 indicate high resistance, whereas negative values indicate lower resistance. We also quantified relative resistance as the difference between total biomass in the replicate and the mean total biomass at 15˚C divided by the mean total biomass at 15˚C (i.e. the loss proportional to the mean biomass at 15˚C).

We did not quantify resilience since in a constant environment temporal variability is in principle closely related to resilience [START_REF] May | Stability and complexity in model ecosystems[END_REF] . To understand potential drivers and mechanisms, we investigated the role of statistical averaging and species asynchrony. Statistical averaging is assessed by the scaling relationship between the species mean biomass and species variance in biomass [START_REF] Tilman | Diversity-stability relationships: statistical inevitability or ecological consequence?[END_REF] . Synchrony was assessed by a metric introduced by Gross et al. (2014) [START_REF] Gross | Species richness and the temporal stability of biomass production: a new analysis of recent biodiversity experiments[END_REF] which calculates the average correlation between the biomass of each species and the total biomass of all other species in the group. We used the R package codyn [START_REF] Hallett | codyn: An r package of community dynamics metrics[END_REF] to calculate asynchrony.

Analyses

Analysis of total biomass used linear mixed effects models [START_REF] Schmid | A guide to analyzing biodiversity experiments[END_REF] (LMMs), with temperature (linear, centred), species richness (log-linear), and their interaction as fixed effects (extended data table 1a). Microcosm identity nested in composition was included as random factor to account for the repeated measurements. Heterogeneity of residuals was accounted for by estimating a composition-specific variance. AIC comparison consistently identified a model with linear temperature effect, as opposed to when temperature was a factor, as most parsimonious (delta AIC of model with nonlinear temperature term was greater than 2 in all comparisons). A simpler analysis of the temporal average of total biomass of each individual microcosm to test the effect of richness, temperature, and their interaction, including a random effect for species composition, yielded quantitatively and qualitatively similar results (see extended data table 1b). Figure 1c shows the estimated effects of LMMs fitted separately for each day with temperature (linear, centred), species richness (log-linear), and their interaction as fixed effects and composition as a random effect. Analyses of temporal stability and resistance did not contain repeated measures and hence only accounted for variation in composition with a random effect. We used the R package nlme for all linear mixed effects models [START_REF] Pinheiro | nlme: Linear and Nonlinear Mixed Effects Models[END_REF] . To calculate associations between resistance and temporal stability among and within richness levels, we calculated reduced major axis regressions (RMA) using the R package lmodel2 [START_REF] Legendre | lmodel2: Model II Regression[END_REF] as both variables potentially contain measurement error. Significance of RMA slopes deviating from zero was assessed by one-tailed permutation tests. All analyses were done in R -the statistical computing environment [START_REF] Core | R: A language and environment for statistical computing[END_REF] .

Review of empirical studies

Based on the review by Donohue et al. (2016) [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF] we obtained a set of studies of resilience, resistance, and temporal variability of ecosystem functions in response to direct or indirect experimental manipulations of diversity. Direct manipulations were defined as changing diversity by adding different sets of species to an experimental plot, jar, or other unit, whereas indirect manipulations induced variation in diversity via the experimental treatment, such as fertilization. We only included studies that performed experimental manipulations. To analyse whether specific mechanisms lead to covariation, we noted the type of mechanism proposed for each of the individual diversity-stability components (extended data table 4). Furthermore, we assessed whether a quantitative or verbal argument was provided (or the mechanisms were not addressed at all) and synthesized the available evidence by vote counting.

Calculating overall ecosystem stability

An interactive web page (https://frankpennekamp.shinyapps.io/Overall_ecosystem_stability_demo/) describes the calculation of ecosystem multifunctionality (also known as overall ecosystem functioning) or overall ecosystem stability and illustrates the following. The calculation requires that values of an ecosystem function (e.g. biomass production) or of a stability component (e.g. resistance to temperature) be converted into a common currency. The threshold approach uses a step mathematical function [START_REF] Byrnes | Investigating the relationship between biodiversity and ecosystem multifunctionality: challenges and solutions[END_REF] ; the averaging approach uses a linear mathematical function (and both equalise relative contributions of different ecosystem functions / stability components) [START_REF] Maestre | Plant species richness and ecosystem multifunctionality in global drylands[END_REF] ; a principal component approach uses a specific linear mathematical function for each ecosystem function or stability component [START_REF] Antiqueira | Warming and top predator loss drive ecosystem multifunctionality[END_REF] ; and Slade et al.

(2017) [START_REF] Slade | The importance of species identity and interactions for multifunctionality depends on how ecosystem functions are valued[END_REF] propose step-like mathematical functions with more or less gradual changes from the lower to higher value. The generalised logistic function (also known as the Richard's function) is flexible enough to give a wide range of shapes of conversion function. If x is the measured variable, and Y is the converted variable, the generalised logistic function is:

Y = A + K -A (C + Qe ,-. ) 0/2
A is the lower asymptote. K is the upper asymptote. B is the gradient. v affects the symmetry, and also the value of y(0). Q affects the value of y(0), i.e. it shifts the function horizontally. C is typically set to 1.

x is a variable, here the value of the measured ecosystem function or stability component. Overall ecosystem stability is then the sum of the standardised and converted stability components OES = f(z(res)) + f(z(ts)), where res is the measured resistance, ts is the measured temporal stability, the function z() subtracts the mean and divides by the standard deviation, and f() is the generalised logistic function. The parameters of f() were A = -1, K = 1, B = 5, v = 1, C = 1 and Q was varied from 10 ,@ to 10 @ . These values were chosen to produce converted stability measures that span the range A to K and to have a relatively threshold-like change from A to K.

Standardisation prior to summation results in overall ecosystem stability with mean of zero, emphasising that the units of valuation here are arbitrary (though generally need not be). Standardisation also implies equal weights for different stability components; weighting of functions needs to be further considered and may be specified according to the specific use cases [START_REF] Gamfeldt | Multiple functions increase the importance of biodiversity for overall ecosystem functioning[END_REF] . Differential weightings, if desired and justified, can be incorporated into the conversions functions. Suggestions regarding the choice of conversion functions for managed systems can be found in Slade et al. 2017 [START_REF] Slade | The importance of species identity and interactions for multifunctionality depends on how ecosystem functions are valued[END_REF] and Manning et al. 2018 [START_REF] Manning | Redefining Ecosystem Multifunctionality[END_REF] .

Unimodal relationships can result from negative covariation among two stability components. How does consideration of more than two components affect the unimodal pattern? While the unimodal relationship is the most pronounced when equal numbers of positive and negative relationships are considered, a unimodal relationship will persist as long as there is at least one opposing stability component (see extended data figure 5).

Analysed data: FP, OP, MP, AT, MS Wrote the first draft: FP, OP Contributed to revisions of the manuscript: all Competing interests: The authors declare no competing interests. Correspondence and requests for materials should be addressed to F.P. or O.P.

Figures

Figure 1: Biomass and stability in experimental microbial communities. Richness increased temporal stability (a1, n=681 independent microcosms), but decreased resistance (a2, n = 567 independent microcosms).

Average (±1 standard error of the mean) biomass for each diversity level (b, n = 12939 microcosm x day combinations). Temporal dynamics of effect sizes (and 95% confidence intervals) of a linear mixed effects model of total biomass showed on average a negative effect of temperature (c1), a positive effect of diversity (c2), and a more negative effect of temperature at higher richness (persistent negative interaction term, c3) (n = 681 independent microcosms per day). Resistance and temporal stability (n = 567 independent microcosms) covaried negatively across richness (d, solid line through centroid of each richness level, n = 6 independent richness levels) but positively within richness levels (dashed lines).

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 1 2 3 4 5 1 2 3 4 5 6

Richness

Temporal stability (mean/STD) Low response diversity was inferred because the biomass of most species decreased or was unaffected by temperature (linear regression between temperature and species biomass; n = 972 species biomass x temperature observations). Consequently, when there were more species, there was greater total biomass, greater temporal stability, but greater biomass loss with temperature increase. Thus, niche complementarity (i.e. effect diversity) likely caused a positive effect of diversity on temporal stability but in the absence of high response diversity, also had a negative effect of diversity on resistance. However, this explanation cannot apply within richness levels, where there was positive covariance among stability components.

a1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -0.

Figure 2 :Figure 3 :

 23 Figure 2: Positive, negative and neutral relationships among resistance, resilience and temporal variability in empirical studies with diversity manipulation. 30 bivariate relationships were reported by 17 independent studies (in addition to this study). Detailed information about individual studies (e.g. code VR-09) is provided in extended data table 3 & 4. Beige regions indicate no covariation. Relative positions within regions are arbitrary and do not indicate relative strengths of relationships. Different colours indicate the effect of diversity on absolute (red) or relative resistance (blue), whereas temporal stability and resilience are shown in black.

1 : 2 :

 12 Richness increased temporal stability across temperatures. a) The stabilizing effect of richness was present across all temperatures, although temperature has a negative effect on mean stability. b) Result table for linear, mixed effects model of log richness, temperature and their interaction on temporal stability supporting the stabilizing effects of richness and the negative effect of temperature on temporal stability (n = 681 independent microcosms). c) Result table for the same analysis as b) but without the monocultures. Results are qualitatively the same, indicating that the relationship between richness and temporal stability is not only driven by the monocultures (n = 580 independent microcosms). Extended data figure The effect of richness on absolute and proportional resistance. a, c) Richness decreased resistance, regardless whether its measured on the absolute or proportional scale. b, d) Result tables of linear, mixed effects model of richness, temperature and their interaction on absolute and proportional richness (n = 567 independent microcosms).

ddata figure 3 :

 3 Extended Niche complementarity and low response diversity were likely caused negative covariance of stability components. Niche complementarity and the resulting increase in total biomass with richness tended to increase temporal stability (figure1in main text). a, b) We found little evidence for an effect of population asynchrony on temporal stability (linear mixed effects model with composition as random effect and log richness and temperature as fixed effects; n = 681 independent microcosms). c, d) In contrast, statistical averaging contributed to stabilization (linear regression between mean species biomass and the variance of species biomass; n = 2077 species mean-variance biomass observations). e)

Acknowledgements

Frederik De Laender and Bernhard Schmid provided valuable feedback on previous drafts of the article. Ian Donohue kindly donated the list of publications from his 2016 review paper. The University of Zurich Research Priority Programme on Global Change and Biodiversity supported this research. Furthermore, funding came from the Swiss National Science Foundation (grant PP00P3_150698 to FA, and 31003A_159498 to OLP).

Code availability

Code to reproduce the analyses and figures is accessible on Github https://github.com/pennekampster/Code_and_data_OverallEcosystemStability (DOI: 10.5281/zenodo.1345557).

Data availability

The experimental data that support the findings of this study are available in Github (https://github.com/pennekampster/Code_and_data_OverallEcosystemStability) with the identifier (DOI: 10.5281/zenodo.1345557).). Source data for figures 1-3 are provided with the paper.

Extended data figure 5:

The effect of aggregating more than two stability components into overall ecosystem stability. The fraction of stability components with negative sign influences whether or not a unimodal pattern will result for a total of 100 stability components. a) A unimodal relationship between diversity and OES will result if at least one stability component is negative. b) However, the strength of the pattern depends on the relative balance of positive and negative relationships. Extended data table 4: Putative mechanisms and type of evidence reported for all bivariate diversity-stability relationships.