Keywords: Powertrain design, hybrid electric vehicles, optimization, sizing, genetic algorithms, dynamic programming, constraint programming, architecture generation. I

Designing a hybrid powertrain remains a complex task. It is an intricate system involving numerous variables that are spread over different levels: architecture, component technologies, sizing, and control. There is currently a lack of frameworks or tools that help in exploring the entire design space and in finding the global optimal solution throughout these levels. This article proposes a systematic methodology that tries to answer a part of this need. Starting from a set of chosen components, the methodology automatically generates all the possible graphs of architectures using constraint-programming techniques. A tailored representation is developed to picture the graphs. They are then transformed into other types of representation (tables describing the connections and the powertrain modes). Based on these representations, the architectures are automatically filtered and the most promising ones are selected. They are automatically assessed and optimized using a specifically developed general hybrid model that calculates the performance and fuel consumption of all the generated architectures. This model is inserted inside a bi-level optimization process: a Genetic Algorithm is used on the sizing and components level, while the Dynamic Programming is used on the control level. A case study is performed and the capability of the methodology is proven.

i. (P)HEV design problem (P)HEV (referring to Hybrid and Plug-in Hybrid Electric Vehicle) powertrains combine a battery electric traction system with a conventional engine-based traction system. Different architectures exist and can be split into 3 main categories: series, parallel and series-parallel [START_REF] Chan | Electric, hybrid, and fuelcell vehicles: Architectures and modeling[END_REF]. Within a category, different variants can be made by changing the location of the powertrain components and by adding gears or gearboxes.

Once the architecture is chosen, different component technologies can be selected, and different sizing can be made. The vehicle operation and fuel consumption will depend on all these choices and on the energy management during the operation.

Figure 1: (P)HEV design problem in its optimization context

In view of this, the design of (P)HEV can be seen as a multiobjective optimization problem that is spread over multiple levels [START_REF] Silvas | Review of Optimization Strategies for System-Level Design in Hybrid Electric Vehicles[END_REF], [START_REF] Van Harselaar | Automated Dynamic Modeling of Arbitrary Hybrid and Electric Drivetrain Topologies[END_REF] (Figure 1). The involved variables can be divided into 3 levels: (1) Architecture level, [START_REF] Chan | Electric, hybrid, and fuelcell vehicles: Architectures and modeling[END_REF] Components technology and sizing level, (3) Control level.

The variables are to be optimized while respecting constraints coming from a system level (as the vehicle performance) and components level (components limitations). The objective functions to optimize can be design related functions (powertrain cost, volume) or design and control related functions (fuel consumption, pollutants emissions).

This design problem is highly challenging due to the high number of involved variables, the numerous conflicting Systematic methodology for architecture generation and design optimization of hybrid powertrains B. Kabalan, E. Vinot, R. Trigui, and C. Dumand T considerations, and the multi-objective, non-linear, non-convex character of the problem. In addition, the assessment of all the combinations of variables (brute force or exhaustive search) is not time feasible.

ii. State of the art

In the literature, this design problem is traditionally tackled by optimizing level 2 (Components technology & sizing) and 3 (Control) without including the level of architecture in the optimization. Few benchmark architectures are arbitrarily selected, manually modeled, optimized on the two levels [START_REF] Chan | Electric, hybrid, and fuelcell vehicles: Architectures and modeling[END_REF] and [START_REF] Silvas | Review of Optimization Strategies for System-Level Design in Hybrid Electric Vehicles[END_REF], and then compared [START_REF] Vinot | Comparison of different power-split hybrid architectures using a global optimization design method[END_REF] [START_REF] Silvas | Review of Optimization Strategies for System-Level Design in Hybrid Electric Vehicles[END_REF]. Some manual topology modifications in the aim of efficiency improvement are also done. In [START_REF] Vinot | Improvement of an EVT-Based HEV Using Dynamic Programming[END_REF] a gearbox is inserted between the Engine (ICE) and the Motor Generator (MG). Different topologies of seriesparallel HEV made by adding gears and gearboxes are studied in [START_REF] Kabalan | Efficiency Improvement of a Series-Parallel Hybrid Electric Powertrain by Topology Modification[END_REF]. The location of the MG can also be changed and improved as in [START_REF] Hofman | Topology Optimization for Hybrid Electric Vehicles With Automated Transmissions[END_REF]. Different power-split topologies can also be considered as in [START_REF] Zhang | Prius(+) and Volt(-): 0% 50% 100% STOP Electric Series Parallel 1MG Parallel 2MG Urban Rural Highway Configuration Analysis of Power-Split Hybrid Vehicles With a Single Planetary Gear[END_REF].

Further, platforms are proposed to compare combinations of vehicle application, architectures, components choice, sizing, and control [START_REF] Gan | Development of a Multi-Architecture and Multi-Application Hybrid Vehicle Design and Management Tool[END_REF] [START_REF] Gan | Multi Architecture Optimization of a Hybrid Electric Vehicle Using Object-Oriented Programming[END_REF]. Sometimes, a general model is used to simulate the few studied architectures. For example, in [START_REF] Chen | Global modeling of different vehicles[END_REF] a general forward Simulink model is developed to assess different architectures by connecting or not some components by means of coefficients equal to a gear ratio or to zero if not connected. General models of power-split architectures are also used in [START_REF] Zhang | Design of Multimode Power-Split Hybrid Vehicles -A Case Study on the Voltec Powertrain System[END_REF] [START_REF] Zhang | Automated Modeling and Mode Screening for Exhaustive Search of Double-Planetary-Gear Power Split Hybrid Powertrains[END_REF] [START_REF] Qin | Simultaneous optimization of topology, control and size for multimode hybrid tracked vehicles[END_REF], see below. In [START_REF] Stroe | Towards a generic control-oriented model for HEV predictive energy management[END_REF], a "matrix" based method with 2 clutches and four possible MGs positions allow to simulate a large class of hybrid vehicles with one gearbox.

In the present paper an original general hybrid model based on the operating modes will be proposed. It will be able to simulate all the hybrid architectures with one engine, two electrical machines and complex or multiple gearboxes.

The architecture level has not been sufficiently explored in the literature. This is because listing and modeling manually all the architectures is infeasible. In addition, it is because of the model complexity and computation time burden associated with the architecture's evaluation and optimization, in case of exhaustive search. Nevertheless, some works started to explore the architecture level by performing automatic generation of architectures instead of sticking to a few known architectures. In [START_REF] Zhang | Design of Multimode Power-Split Hybrid Vehicles -A Case Study on the Voltec Powertrain System[END_REF][15] [START_REF] Zhang | Automated Modeling and Mode Screening for Exhaustive Search of Double-Planetary-Gear Power Split Hybrid Powertrains[END_REF], a generation and comparison of power-split architectures is performed based on system matrix generation and analysis. The principle is to use one or two planetary gears (PG) linked to the other components by clutches which can be closed or open. It allows the assessment of different topologies and operating modes depending on the state of the clutches. However, only a brute force comparison is performed with the components and gear-ratio of the Volt-2 and only some downsizing of the MGs is presented. In [START_REF] Qin | Simultaneous optimization of topology, control and size for multimode hybrid tracked vehicles[END_REF] a new model of tracked vehicle with three PG is developed and optimized sizing is performed using genetic algorithm (GA) and near optimal control method. In most of these studies, the efficiency of the gears seems not to be taken into account, thus energetic comparison could not be completely fair.

In [START_REF] Silvas | Functional and Cost-Based Automatic Generator for Hybrid Vehicles Topologies[END_REF], the automatic generation of graphs of hybrid architectures is proposed using constraint-programming, without an evaluation of the architectures. Automated physical modeling and filtering is then added in [START_REF] Wijkniet | Modified Computational Design Synthesis Using Simulation-Based Evaluation and Constraint Consistency for Vehicle Powertrain Systems[END_REF] without component technology & sizing optimization, and while neglecting the transmission losses. These works will be recalled later and compared to the proposed methodology.

In this paper, the authors propose a complete methodology allowing the optimal design of HEV architectures from topology level to the optimal sizing and control. Using a certain amount of components, all the possible topologies are generated, filtered and optimally sized (components and gears). An original model allows assessing the fuel consumption while taking into account the efficiency of all the components and without the need todevelop a new model when changing the involved components.

iii. Proposed methodology Figure 2: Overview of the proposed methodology

The proposed methodology (Figure 2) starts by a definition of the design inputs (A part in Figure 2): the components that can be used, and the specifications (vehicle requirements, objective functions).

Then an automatic generation of all the possible architectures (B part in Figure 2) is done using constraint-programming. This method is similar to the one used in [START_REF] Silvas | Functional and Cost-Based Automatic Generator for Hybrid Vehicles Topologies[END_REF]; however, in this work the problem is differently defined, more conveniently implemented using Python environment, and easily solved (chapter III). In this step, the architectures are represented as graphs using a novel graphical representation (chapter II) suitable for new-trend hybrid gearboxes that cannot be generated and studied using the methodologies in the previously cited works [START_REF] Zhang | Design of Multimode Power-Split Hybrid Vehicles -A Case Study on the Voltec Powertrain System[END_REF][17] [18][19]. Thousands of generated graphs are commonly expected. Some filtering steps could be added in order to select tens of graphs which correspond to the most promising architectures. This filtering (chapter IV) is done based on novel representations developed in this work.

The most promising architectures should be evaluated and optimized. A General Hybrid Model (C part in Figure 2) is developed for this purpose. It is inserted inside a bi-level optimization technique (D part in Figure 2) on the levels of control and component technology & sizing. To solve the linked problem of sizing and control (level 2 and 3, figure 1) the most used method seems to be a joined optimization of component sizes and control. This is often performed using an outer loop of optimization based on stochastic algorithm and an inner control loop which can be based on optimal control as dynamic programming (DP) method or on near optimal control method [START_REF] Qin | Simultaneous optimization of topology, control and size for multimode hybrid tracked vehicles[END_REF]. For example, in [START_REF] Silvas | Comparison of Bi-level Optimization Frameworks for Sizing and Control of a Hybrid Electric Vehicle[END_REF] and [START_REF] Zou | Combined optimal sizing and control for a hybrid tracked vehicle[END_REF] GA and DP are used. In [START_REF] Pourabdollah | Optimal sizing of a series PHEV: Comparison between convex optimization and particle swarm optimization[END_REF] particle swarm optimization are tried with DP. Other techniques are also tried, using Pontryagin's minimum principle for sizing and control [START_REF] Song | Component sizing optimization of plug-in hybrid electric vehicles with the hybrid energy storage system[END_REF] or even gradient based technique for the outer loop [START_REF] Poline | Simultaneous optimization of sizing and energy management-Application to hybrid train[END_REF]. However, due to the nonconvexity of the problem, GA associated to DP seems to be the most robust method to tackle the sizing and control problem [START_REF] Silvas | Comparison of Bi-level Optimization Frameworks for Sizing and Control of a Hybrid Electric Vehicle[END_REF], and performs a fair comparison while guaranteeing the convergence and the optimal control.

The developed General Hybrid Model is inserted inside a bilevel optimization technique (GA and DP) and is able to automatically assess and optimize all the generated architectures, an achievement not yet seen in the literature. The architectures are finally compared based on the chosen objective functions.

iv. Contribution and outline

The contributions of this paper are: 1-Exploration of the architecture level instead of selection of few known architectures for comparison. This is done thanks to the automatic generation of architectures. 2-Development of a general hybrid model that can allow the assessment of all the generated graphs and therefore can connect the exploration work on the architecture level with the optimization work on the components technology & sizing and control levels. 3-Combining the exploration of the architecture level with the automatic multi-objective optimization. This was not previously achieved in the literature.

In the following, the representation used to visualize the architectures is revealed in chapter II. The automatic generation of architectures is described in chapter III and the automatic filtering in IV. The automatic optimization done using the general hybrid model is explained in chapterV. Finally, a case study is performed in VI to prove the capability of the methodology before concluding.

II. PROPOSED REPRESENTATION

The methodology proposes to generate all the possible architectures that can be made when connecting a certain number of defined powertrain components and defined connecting elements. For that, these components and generated architectures need to be visualized.

In the literature, two remarkable representations are found in [START_REF] Silvas | Functional and Cost-Based Automatic Generator for Hybrid Vehicles Topologies[END_REF], and [START_REF] Masfaraud | Automatized gearbox architecture design exploration by exhaustive graph generation[END_REF]. Inspired by these two works, a graphical representation is proposed to be used in this article (Figure 3).

i. Overview of the graphical representation

A powertrain is visualized as a graph. This graph is made of nodes and connections. The nodes are of two types (bipartite graph) [START_REF] Masfaraud | Automatized gearbox architecture design exploration by exhaustive graph generation[END_REF]: components (circles in colors) and connectors -i.e. shafts -(small circles in black). No direct connection can be done between the components. Connections are only done through the shafts.

The components nodes are:  Powertrain components: the Internal Combustion Engine (ICE), the Motor Generators (MGs), the Final Drive output (FD)  Gearbox components: the clutches and the 'synchro' unit (Synchronizer + 2 gear pairs) ii. The 'synchro' unit in details More details are given here on the synchro unit concept because it gives the representation its uniqueness and appropriate level of details. Such units are seen in manual gearboxes and consist of a synchronizer and 2 pairs of gears (Figure 4). In each pair of gears, one gear (the green) is fixed to a shaft and one (the blue) is free spinning and connects to another shaft through the synchronizer. The synchronizer can be in 3 states: (1) in the center and not connected to any freespinning gear: position Neutral, (2) to the left and connected to the left gear: connection with gear ratio 1, (3) to the right and connected to the right gear: connection with gear ratio 2.

Figure 4: The synchro unit

This 'synchro' unit is visualized as an entity of 3 red nodes (in Figure 4: S_0: central node, S_R: right node, S_L: left node) and has three external connections. One connection is to the input shaft and the two other connections can be to a same output shaft (case of conventional manual gearboxes) or to 2 output shafts (case of dedicated hybrid gearboxes newly seen, as in [START_REF] Fremau | Transmission pour vehicule automobile a propulsion hybride et procédé de commande associé[END_REF][27] [28][29]).

This gives the proposed representation a level of details appropriate to generate these new-trend dedicated hybrid gearboxes that cannot be generated by [START_REF] Silvas | Functional and Cost-Based Automatic Generator for Hybrid Vehicles Topologies[END_REF] without making the problem too complex like in [START_REF] Masfaraud | Automatized gearbox architecture design exploration by exhaustive graph generation[END_REF]. In fact, in [START_REF] Silvas | Functional and Cost-Based Automatic Generator for Hybrid Vehicles Topologies[END_REF] the gearbox is considered as one node that can have two connections (1 input, 1 output) and that can be used once (as in conventional powertrains). There is no possibility to innovate inside the gearbox. In [START_REF] Masfaraud | Automatized gearbox architecture design exploration by exhaustive graph generation[END_REF], all the elementary gearbox components are visualized (the synchronizer and every single gear). This gives a great degree of freedom inside the gearbox, but makes the problem more complex if one is only interested in the energetic aspect of the gearbox and is not in its volumetric and compactness aspect.

iii. Example of architectures

Below are examples on architectures visualized using the proposed representation:  A simple series-parallel architecture with gears for the electric machines: It should be mentioned that in the proposed representation no planetary gear component was used at the time, but it can be easily added in the future. The general model (section V.i) will be adapted.

III. AUTOMATIC GENERATION

In this section, the automatic generation of architectures is explained. This involves solving a Constraint Satisfaction Problem (CSP) presented hereafter.

i. Problem definition 1) Variables and domains:

An architecture is visualized as a graph: a set of nodes V connected by a set of edges E. The graph can be described using an adjacency matrix: 𝑛 by 𝑛 matrix where the columns and the rows are the node 𝑉 0 , … 𝑉 𝑛-1 , and each cell at index (𝑖, 𝑗) takes a 0 or 1 value describing the absence (0) or presence [START_REF]Europe accelerates the transition to clean mobility: Co-legislators agree on strong rules for the modernisation of the mobility sector[END_REF] of connection between the node 𝑉 𝑖 and 𝑉 𝑗 .

In this work, the powertrain is considered to include a maximum of: 1 ICE, 2 MGs, 1 FD, 4 shafts, 3 clutches, and 3 synchro units. With these components quite all new compact parallel or series-parallel architectures for passenger vehicles can be generated, such as presented in [START_REF] Fremau | Transmission pour vehicule automobile a propulsion hybride et procédé de commande associé[END_REF] and [START_REF] Fremau | Hybrid transmission with offset electric machine and method for controlling gear changes[END_REF]. Topologies with multiple gearboxes can be generated including architectures where the MGs are connected to some internal gears of the gearbox and not just to the input or the output shafts of the gearbox. Architectures with more components can be encountered in heavy vehicles but are not relevant (cost issue) for passenger vehicles. The method can be adapted to more components but will be more time consuming.

The corresponding adjacency matrix is shown in Figure 7. In blue and green are the powertrain components' nodes (ICE, MGs, FD), in red are the mechanical transmission components' nodes (clutches, synchronizers), and in grey the shafts nodes (X1 to X4).

The matrix has 23 × 23 = 529 cells that can take 2 values (0 or 1). The values of the variables are related between themselves through some constraints. The involved CSP can be defined as follows:

 Number of variables: 529  Domain of all variables: {0, 1}  Constraints: see next section  Target: get all the feasible solutions The number of solutions without constraints is:

2 529 = 1.7574 × 10 159 (1)

2) Constraints:

Any solution should respect some constraints regarding the components connection, mechanical feasibility, powertrain functionality and the non-repeatability of the components. A summary of the used constraints is shown in Table 1.

Table 1: Summary of the problem constraints

These constraints reduce the design space. Before adding these constraints, the variables are 529 (violet rectangle in Figure 7) and have unknown values. The number of solutions is 2 529 .

The C000 constraint (of Table 1) imposes the matrix to be symmetric because the graph of mechanical connections is an undirected graph: node 1 is connected to node 2 is equivalent to node 2 is connected to node 1. C001 imposes the nodes to be connected to the shaft nodes only, enforcing the variables of other connections to be 0 (except some values kept 1 to visualize the mandatory connection between the nodes of the synchro unit). This is seen in Figure 7 visualizing how the C000 and C001 are affecting the adjacency matrix. The variables with unknown value are reduced from 529 to 76 and the number of solutions from 2 529 to 2 76 . The other constraints are also reducing the design space. They are functionality and nonrepeatability constraints (for example: in C008 it is useless to have 2 clutches in parallel). Another case of constraint is C003 which imposes that, in order to forbid direct connection between ICE and FD, the value at index (0, j) and the value at index (3, j) cannot be both equal to 1, for a column index j between 19 and 22. The problem needs to be implemented in an environment and solved. Constraint programming techniques will be used. Python-Constraint module [START_REF] Niemeyer | Python-Constraint[END_REF] on Python environment was selected. It offers backtracking capabilities and forward checking capabilities (constraint propagation). It succeeded in solving our problem in a relatively fast way. Figure 8 presents an example with a reduced possible number of components.

The output solutions are represented by their adjacency matrices. They need to be visualized according to our proposed representation. For this, the python NetworkX [START_REF] Hagberg | Exploring network structure, dynamics, and function using NetworkX[END_REF] module was used. Between the 132 generated graphs in this example, some are unique but some are redundant. Two examples of redundant graphs are shown in Figure 9. The 2 graphs are the same except a detected interchangeability between the nodes of the synchro unit (S1_R and S1_L).

All the redundant graphs but one need to be removed. This is done in the automatic filtering steps.

IV. AUTOMATIC FILTERING

Even after considering the previously presented constraints (Table 1), the number of possible solutions remains very high and all the generated architectures are not necessarily efficient or different. Therefore, different steps of filtering are needed. They are explained in this chapter. Some graphs were found to be redundant or isomorph due to a symmetry problem: orientation of the synchro unit, interchangeability between the nodes of the synchro unit and interchangeability between the electric machines. The 2 graphs in Figure 9 are examples of the interchangeability between the nodes of the synchro unit (S1_R and S1_L). The 2 graphs have the same energetic model leading to same energy efficiency. There is no need to consider both graphs in the methodology and hence one of them is to be removed.

A new representation called '0ABC Table ' was developed to detect these redundancies. It is a table describing the type of 1stage connection between all the powertrain components. A 1stage connection means that a maximum of 2 shafts exists between the 2 considered components. Figure 10 shows the '0ABC Table' of the redundant graphs of Figure 9. The 1 st line of the 0ABC Table indicates that the ICE can be connected to FD through a synchro. The 2 nd line: between MG1 and FD, 1 connection through synchro is possible. The 3 rd line describes the connection between MG2 and FD, but it is set to 0 since there is no MG2 in this architecture. 4 th line: 2 connections through synchro can be done between ICE and MG1. 5 th and 6 th lines are 0 since MG2 is not present.

Graphs having the same 0ABC table are redundant graphs. The decision is to keep the 1 st graph and remove the others.

Figure 10: 0ABC table of the redundant graphs in Figure 9 C003: Values at a same column index j cannot be both equal to1

ii. Filtering Step 2: Modes Table

Each of the kept non-redundant graphs has different states caused by the different combination of states of the actuators. A graph equipped with 𝑁_𝑐𝑙𝑢𝑡𝑐ℎ𝑒𝑠 and 𝑁_𝑠𝑦𝑛𝑐ℎ𝑟𝑜𝑠 will have 2 𝑁_𝑐𝑙𝑢𝑡𝑐ℎ𝑒𝑠 × 3 𝑁_𝑠𝑦𝑛𝑐ℎ𝑟𝑜𝑠 states. The developed tool on Python generates these states graphs and automatically detects the corresponding powertrain mode, based on the components present and which components are connected. Nine modes are possible:

(1) ICE only: ICE alone is connected to the wheels (2) Pure 1MG: electric mode, only 1 MG connected to the wheels (3) Pure 2MG: electric mode, 2 MGs connected to the wheels (4) Parallel 1MG: ICE and 1 MG connected to the wheels (5) Parallel 2MG: ICE and 2 MGs connected to the wheels (6) Series: 1 MG is connected to the wheels and the ICE is connected to another MG (7) Neutral: the wheels are connected to none of the components, and no connection between the other powertrain components (8) Stand-still charging 1MG: the wheels are connected to none of the components; the engine is however connected to 1 MG (9) Stand-still charging 2MG: the wheels are connected to none of the components; the engine is however connected to 2 MGs These modes were imagined to cover the possibilities of energy paths. The option of 'undetected mode' was added in case a state graph did not fit to the 9 possible modes. However, this case was never encountered proving that the 9 possible modes do cover all the possible modes.

For each graph, a 'Modes Table ' is created and lists the number of each of the 9 possible modes available in the graph (Figure 11).

Figure 11: Modes Table

At this level, the graphs can be compared based on their Modes Table . In addition, constraints can be added on the minimum and maximum number of modes required. An example is shown in Figure 11. The Modes Table of the considered architecture of the example has 0 "ICE only" mode. This could be used for example as a constraint and will dismiss this architecture as the minimum number of "ICE only" modes required in this example is 1.

iii. Toward assessment: Modes Table +

More information is needed to assess the architectures and is linked to the efficiency of the powertrain in each of the listed modes. For this reason, a Modes Table + is created (Table 2) and will be used in the assessment (part V.i). For each available mode, a description of the paths between all the components is done by calculating:

 The global efficiency of the path  The global gear ratio of the path Each node in the graph was already assigned an attribute ratio and efficiency. The global gear ratio or efficiency of a path is calculated as a multiplication or division of the attributes of the nodes that are on the path. An example of how some modes are described in Modes Table + is shown in Table 2.

For instance, the ICE only mode (first line of Table 2) is described using 2 values: K_ICE_FD and η_ICE_FD. K_ICE_FD is the global ratio of the ICE to FD path and is calculated as a multiplication of all the ratios of the nodes that exist on this path. η_ICE_FD is the global efficiency of the ICE to FD path and is also calculated as a multiplication of the efficiencies of the nodes on the path. For the other modes, all the paths between the involved components are described in a similar way. When only 1 MG is involved in a mode (example Pure 1MG), a variable MG is given a value of 1 or 2 to define which MG is connected to the wheels. The Modes Table + include the information needed for the energetic assessment of the architectures. This automatic assessment and optimization of the generated architectures will be the content of the next chapter.

Table 2: Example of modes description in Modes Table +

V. AUTOMATIC OPTIMIZATION

The generated architectures need to be assessed, optimized and compared. Therefore, energetic models that calculate the fuel consumption and the performance of the powertrains are required. When the number of assessed architectures is limited, the traditional way is to manually develop one model per architecture. However, when the automatic generation of architectures is included in the methodology, manually modeling all these generated architectures becomes infeasible.

To solve this, automatic modeling [START_REF] Goos | Automatic Evaluation and Optimization of Generic Hybrid Vehicle Topologies using Dynamic Programming[END_REF] or generic transmission models [START_REF] Van Harselaar | Automated Dynamic Modeling of Arbitrary Hybrid and Electric Drivetrain Topologies[END_REF] can be used. In this work, a specific general model for (P)HEV is developed.

i. General hybrid model

This model includes the 9 possible modes that the generated architectures can have. A visualization of the model is shown in Figure 12. No real transmission exists in this model but it can represent all the generated architectures. Even if the energy path is not real, it is perfectly equivalent to the "real" generated topologies and allows assessing the vehicle performance or the fuel consumption in an exact manner. All the components are connected by simplified virtual links that have 1 gear with a ratio and efficiency. These are the global ratio and efficiency of the considered power path. This information is found in the Modes Table + of the architecture.

Figure 12: General Hybrid Model

Depending on the selected mode, the energy path is activated or not. Examples of 2 modes are shown in Figure 13. Each of the possible modes in the architecture is a function called using the description of the mode found in Modes Table+.

Figure 13: Example of modes inside the general hybrid model

The general hybrid model includes a function for the dynamic performance calculation of the vehicle and a function for the fuel consumption calculation using Dynamic Programming (DP), see part ii. In each of these 2 scripts, the calculation requires an instantaneous calling of all the available modes. All the available modes in Modes Table + are called with the specific information that is also found in Modes Table +.

ii. Bi-level optimization process

The general hybrid model has components & sizing variables that need to be decided when the powertrain is being designed, and control variables that need to be decided instantaneously when the powertrain is being operated. To optimize these variables, the model is inserted inside a bi-level optimization technique following the nested approach [START_REF] Silvas | Review of Optimization Strategies for System-Level Design in Hybrid Electric Vehicles[END_REF]: the design (components technology & sizing) in an upper level and the control in a lower level (Figure 14).

1) Upper level design optimization:

A Genetic Algorithm (GA) manages the design level and makes the choice of design candidates to assess (set of values for the design variables). Each design candidate will pass a set of performance tests (for example: acceleration 0 to 100 km/h, acceleration 80 to 120 km/h, and maximum speed). For that, the performance function is called. Inside this function, the vehicle starts from stationary position and at each time step we search for the torque that can be delivered to the wheels. For this, all the available modes are called using their specific information found in Modes Table+. The mode that can deliver the maximum torque on the wheels is selected.

Figure 14: Bi-level optimization process and function calling

For the candidates that pass the performance test, a calculation of the objectives is done. GA iterates this process and performs an exploration of this design spaceusing an evolution process (Selection, Crossover, and Mutation) to produce better generations. This process yields to a Pareto front presenting the tradeoff between the objective functions.

2

) Lower level control optimization:

A calculation of the objectives is needed on the upper level of design optimization. When the objective function is only design related (powertrain cost, volume, sizing), no vehicle operation is needed and no control optimization is needed. However, in the case of an objective like the fuel consumption of the vehicle, the vehicle should be simulated and its control is to be optimized.

The vehicle is considered to operate on a driving cycle and Dynamic Programming DP [START_REF] Bellman | Dynamic Programming[END_REF][34] solves the control problem. This problem can be represented as a graph of battery State Of Charge SOC versus time (the SOC being the only state variable in the problem). The objective is to go from an initial SOC to a same final SOC (charge sustaining) while choosing the control that ensures the minimum cumulative fuel consumption at the end of the driving cycle. DP solves this problem in 3 steps: graph construction (space limit of the problem), edge cost calculation and finally graph solving. For the first 2 steps, the limits function and the edge cost function of the general hybrid model are called. The final SOC constraint is ensured by the construction of the searching space. The instantaneous fuel consumption cost of the edges is assessed by means of the general model presented (parti i), and DP finds the optimal control to minimize the global fuel consumption on the cycle. For more information on DP applied to HEVs see for example [START_REF] Vinot | Improvement of an EVT-Based HEV Using Dynamic Programming[END_REF]. Inside these functions, all the available modes are to be called instantaneously. Based on Modes Table +, all the available modes are called using their specific information. Unlike conventional ways, DP is operating here with no previous knowledge of the architecture.

The GA, DP and all the scripts are fixed. Changing Modes Table + and repeating the process is enough to simulate and optimize all the generated architectures.

VI. CASE STUDY

In this chapter, an application of the entire methodology is performed to show its capability. Starting from the powertrain components of a Series-Parallel powertrain, all the feasible architectures are generated and filtered. The most promising are automatically optimized and compared.

i. Input components and generated graphs

The starting components that are injected in the automatic generation tool are the components of a series-parallel architecture that has a 2-speed gearbox (Figure 15): 1 ICE, 1 FD, 2 MGs, 4 shafts (black dots), 2 clutches (red dots) and 1 synchro unit. We will be interested in architectures having at least 1 ICE only, 1 electric, 1 series and 1 parallel mode. In addition, the architectures should allow battery charging at standstill (cf.VI.ii).

The tool generated 480 graphs of architectures in 807.1 seconds. The eleven architectures finally selected (see below) are presented in Figure 16.

ii.Automatic filtering and selection of architectures

The 480 generated graphs are filtered in 3 steps: At the end of the filtering steps, 11 graphs of architectures remained (Figure 16). These are the most promising architectures that are optimized and compared in the next section.

iii. Automatic assessment and comparison

The Modes Tables + are passed to the general hybrid model.

1) Simulation parameters

The components used correspond to a middle class HEV [START_REF] Kabalan | Efficiency Improvement of a Series-Parallel Hybrid Electric Powertrain by Topology Modification[END_REF]. The components are sized through a scaling technique [START_REF] Kabalan | Efficiency Improvement of a Series-Parallel Hybrid Electric Powertrain by Topology Modification[END_REF] [START_REF] Vinot | Improvement of an EVT-based HEV using dynamic programming[END_REF]. Their initial main characteristics are presented Table 3. The engine specific fuel consumption and MGs efficiencies maps can be found in [START_REF] Kabalan | Efficiency Improvement of a Series-Parallel Hybrid Electric Powertrain by Topology Modification[END_REF]. The sizing variables are scaling factors applied to the maximum power of the ICE, MG1 and MG2, the modules number in the battery and all the involved gear ratios [START_REF] Vinot | Improvement of an EVT-Based HEV Using Dynamic Programming[END_REF]. The scaling factor also impacts the weight of the components (with the assumption of weight proportional to maximum torque). The initial weights are presented in Table 3. 𝛼, 𝛽 and 𝛾 are coefficients calculated from the mean traveled distance by the French population in urban, rural road and highway conditions. These values are respectively 0.4, 0.3 and 0.3 [START_REF]Eurostat statistic on population displacement[END_REF].

The fuel consumption will be a good indication of the CO2 emission. Further works can also consider other pollutants emissions (CO, NOx, HC). However it will complicate the control problem if a tradeoff between CO2 and pollutant is also considered as for example in [START_REF] Guille Des Buttes | Energy management strategy to reduce pollutant emissions during the catalyst light-off of parallel hybrid vehicles[END_REF]. The battery number is a good indicator of the onboard electrical power and energy and of the cost of the hybridization. A global cost can also be considered [START_REF] Kabalan | Optimal Design and Sizing of Through-The-Road Hybrid Vehicle Powertrain[END_REF][15] even if it is difficult to evaluate with accuracy.

The multi-objectives NSGA-II algorithm will then provide the Pareto front of the optimization problem formulated as follows:

{ min(𝐹𝐶𝑜𝑛𝑠) min(𝑁𝑏𝑎𝑡𝑡) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑡ℎ𝑒 𝑓𝑜𝑙𝑜𝑤𝑖𝑛𝑔 𝑐𝑜𝑛𝑠𝑡𝑟𝑎𝑖𝑛𝑡 C1:𝑡0->100< 10.1 s; C2:𝑡80->120< 7.5 s; C3:Vmax> 180 km/h FCons is computed as a weighted average of the fuel consumption in charge sustaining mode in urban, rural road and highway conditions. The ARTEMIS European driving cycles are used to simulate these three conditions and the control is optimized on each of them using DP. One DP problem is solved (part V.ii.2) for each cycle of FCons.

2) Results

The 11 architectures are optimized. The Pareto fronts representing the performance of each architecture according to the 2 objectives are shown in Figure 17. The Pareto fronts shown in Figure 17 prove that for a same number of battery modules, the architectures can have a 10% to 14% difference in fuel consumption. For a same fuel consumption, the architectures can have 40% to 60% difference in the number of battery modules in series. 4-Another output of the methodology is the possibility of examining closely the optimal sizing and the powertrain operation (energy management, components operation, modes choice…) for all the architectures on all the simulated driving cycles. An example is presented hereafter. The Pareto points corresponding to 28 Number of battery modules for architectures 9 (one of the best), 1 (intermediate) and 2 (worse) are considered in this example. The sizing is presented in Table 4, along with the fuel consumption in each of the driving conditions. Architecture 9 has the least fuel consumption in all driving conditions. It has also downsized MG1 and ICE compared to architecture 2. However, architecture 1 presents the most downsized electric powertrain components (MG1 and MG2), with a slightly upsized engine.

The operation of architecture 9 on the rural driving cycle is examined in Figure 18 and Figure 19. The charge sustaining operation is evident and frequent engine power cut is visible. The percentage of time choice of each mode in each driving condition is presented in Figure 20. The electric driving has the major shares in urban and rural, with less occurrence in highway. The parallel 1MG modes usage increases from urban to rural and highway. The parallel 2MG modes are absent in urban and minorly used in rural and highway. The series mode is also minorly used in all the driving conditions. The importance of having a series mode is to be questioned in later analysis.

VII. CONCLUSION

The paper proposes a systematic methodology for the generation and design of hybrid powertrains. The representation used and the different parts of the methodology were explained: the automatic generation of architectures, the automatic filtering and the automatic optimization. A case study was also done and proved the capability of the methodology. The main contributions of this work are (1) the developed tool that can generate and filter hybrid architectures based on the manufacturer constraints and requirements, and (2) the general hybrid model that can simulate and optimize all the generated architectures. The main limitations are the exclusion of powersplit architectures and the optimization time for 1 architecture (8 to 10 hours). These issues will be addressed in the future. Improving the codes has already started and could reduce the optimization time by half. The component models, sizing technique and choice will also be reviewed and improved. Other objectives will be also considered (compactness, pollutants, battery aging...). An example of comprehensive analysis of the results (optimal sizing, powertrain operation, modes choice) has been given in this article, but more deep analysis will be carried out in the future in order to assess in detail the effectiveness of the best proposed solutions.

Figure 3 :

 3 Figure 3: The proposed representation

Figure 5 : Example 1

 51 Figure 5: Example 1 Gears are not shown in our representation. However, each component edge and each of the right and left nodes of the 'synchro' unit could have an attribute gear ratio.  A parallel architecture with a 2-speed gearbox:

Figure 6 :

 6 Figure 6: Example 2

Figure 7 :

 7 Figure 7: The adjacency matrix after the addition of constraints ii. Problem implementation and solvingThe problem needs to be implemented in an environment and solved. Constraint programming techniques will be used. Python-Constraint module[START_REF] Niemeyer | Python-Constraint[END_REF] on Python environment was selected. It offers backtracking capabilities and forward checking capabilities (constraint propagation). It succeeded in solving our problem in a relatively fast way. Figure8presents an example with a reduced possible number of components.The output solutions are represented by their adjacency matrices. They need to be visualized according to our proposed representation. For this, the python NetworkX[START_REF] Hagberg | Exploring network structure, dynamics, and function using NetworkX[END_REF] module was used. Between the 132 generated graphs in this example, some are unique but some are redundant. Two examples of redundant graphs are shown in Figure9. The 2 graphs are the same except a detected interchangeability between the nodes of the synchro unit (S1_R and S1_L).All the redundant graphs but one need to be removed. This is done in the automatic filtering steps.

Figure 8 :Figure 9 :

 89 Figure 8: Example of solving

  0ABC filtering: Removed 417 redundant graphs and kept 63 non-redundant graphs. Modes filtering: Architectures are required to have: -1 ICE only mode -1 (Pure 1MG or 2MG) mode -1 Series mode -1 (Parallel 1MG or 2MG) mode -1 (Series or Stand-still charging 1MG or Stand-still charging 2MG) mode This step dismissed 46 graphs and 17 remained.

Figure 15 :Figure 16 :

 1516 Figure 15: Input components and solving

  The battery is modeled by a classical equivalent electric serial model: open circuit voltage in series with a resistance. The resistance value changes depending on charge or discharge operations. Two objective functions are considered: (1) Number of battery modules in series (Nbatt), (2) Charge sustaining fuel consumption: 𝐹𝐶𝑜𝑛𝑠 = 𝛼. 𝐹𝐶𝑜𝑛𝑠 𝑢𝑟𝑏𝑎𝑛 + 𝛽. 𝐹𝐶𝑜𝑛𝑠 𝑟𝑢𝑟𝑎𝑙 + 𝛾. 𝐹𝐶𝑜𝑛𝑠 ℎ𝑖𝑔ℎ𝑤𝑎𝑦 (1)

Figure 17 :

 17 Figure 17: Pareto fronts at the end of the optimization Those results prove that the methodology succeeded in: 1-Automatically generating all the graphs of possible architectures. 2-Automatically filtering these architectures and selecting some of them for optimization and assessment. 3-Automatically optimizing these architectures and comparing them based on the 2 chosen objective functions.The Pareto fronts shown in Figure17prove that for a same number of battery modules, the architectures can have a 10% to 14% difference in fuel consumption. For a same fuel consumption, the architectures can have 40% to 60% difference in the number of battery modules in series. 4-Another output of the methodology is the possibility of examining closely the optimal sizing and the powertrain operation (energy management, components operation, modes choice…) for all the architectures on all the simulated driving cycles. An example is presented hereafter.

Figure 18 :Figure 19 :

 1819 Figure 18: SOC variation in rural driving

Figure 20 :

 20 Figure 20: Mode choice in percentage of time in urban, rural and highway driving

Table

Table 3 :

 3 Main initial characteristics of the components

	Vehicle Weight		1360 kg
	Engine	weight	100 kg
		max power	57kW @5000rpm
	MG1 and MG2	weight	55 kg
		max power	50 kW
		max torque	400Nm @1200 rpm
	Battery	Module weight	1.75 kg
		Module energy	46.7 Wh

Table 4 :

 4 The sizing of the Pareto points for 28 battery modules

ACKNOWLEDGMENT

The authors acknowledge the contribution of Yuan CHENG (chengyuan.hit@gmail.com) to this work when he was with Groupe PSA.


Vincent Baillard
email: vincent.baillard@univ-lorraine.fr

Marie Laure Delignette-Muller

Cécile Sulmon

Anne-Kristel Bittebiere

Cendrine Mony

Ivan Couee

Gwenola Gouesbet

Simon Devin

Elise Billoir

How does interspecific competition modify the response of grass plants against herbicide treatment? A hierarchical concentration-response approach

Keywords: hierarchical modelling, environmental risk assessment, isoproturon, herbicide, biotic interaction, multi-stress 1

Ecological interactions are rarely taken into account in environmental risk assessment. The objective of this work was to assess how interspecific competition affects the way plant species react to herbicides and more specifically how it modifies the concentration-response curves that can be built using ecotoxicological bioassays. To do this, we relied on the results of ecotoxicological bioassays on six herbaceous species exposed to isoproturon under two conditions: in presence and in absence of a competitor. At the end of the experiments, eleven endpoints were measured. We modelled these data using a hierarchical modelling framework designed to assess the effects of competition on each of the four parameters of the concentration response curves (e.g. the level of response at the control or the concentration at the inflexion point of the curve) simultaneously for the six species. The modelled effects could be of three types, 1) competition had no effect on the parameter, 2) competition had the same effect on the parameter for all species and 3) competition had a different effect on the parameter for each species. Our main hypothesis was that different species would react differently to competition. Results showed that about a half of the estimated parameters showed a modification under competition pressure among which only a fourth showed a species-specific effect, the three other fourth showing the same effect between the different species. Our initial hypothesis was thus not supported as species

Introduction

Environmental Risk Assessment (ERA) is mainly based on data from monospecific bioassays. The obtained effect concentrations are usually divided by an Assessment Factor (AF) to compensate for processes not taken into account [START_REF] Aagaard | Guidance on tiered risk assessment for plant protection products for aquatic organisms in edge-of-field surface waters[END_REF] such as biotic interactions, intraspecific variability or interlaboratory variability. The value of the applied AF depends on the method applied to calculate predicted no effect concentrations, the quality of the used data, the ecological relevance of these data and the type of ecosystem considered. For ERA in freshwater, the AF value is for example comprised between 1 (for data with high ecological relevance and method including numerous species and modelling process such as species sensitivity distributions) and 1,000 (for basic procedure keeping the lowest critical effect concentration value from three bioassays) [Amiard and Amiard-Triquet].

Interspecific competition is one of the biotic interactions that can have an important effect on plant development. It is defined by Aschehoug et al. as "the ability of individuals to usurp resources or otherwise suppress their neighbor's fitness and include both resource and interference competition" [START_REF] Aschehoug | The Mechanisms and Consequences of Interspecific Competition Among Plants[END_REF]. Even if this biotic constraint has been largely studied [START_REF] Aschehoug | The Mechanisms and Consequences of Interspecific Competition Among Plants[END_REF], studies of its interactions with chemical stress remain scarce. For example, Boutin et al. studied how sublethal doses of two herbicides (glyphosate and metsulfuron methyl) could modify the competitive interactions between two plants (Centaurea cyanus L. and Silene noctiflora L.) [START_REF] Boutin | Effects of Sublethal Doses of Herbicides on the Competitive Interactions Between 2 Nontarget Plants, Centaurea cyanus L. and Silene noctiflora L[END_REF]. The results showed that competitive ability of S. noctiflora was negatively affected by glyphosate, leading to a competitive advantage of C. cyanus.

Concerning metsulfuron methyl, interspecific competition was found to have an important effect on the herbicide effects on seed production. Another study, from Damgaard et al. tested the effects of glyphosate on growth and competitive effect of two grass species (Festuca ovina and Agrostis capillaris) in a semi-field experiment [START_REF] Damgaard | The effect of glyphosate on the growth and competitive effect of perennial grass species in semi-natural grasslands[END_REF]. The results showed that competitive relationships were modified as the competitive effect of F. ovina increased whereas that of A. capillaris decreased with increasing doses of glyphosate. It is important to emphasise that the study of the interaction between interspecific competition and chemical stress encompasses two questions: whether the competitive ability is affected by the chemical stress or whether the tolerance to the chemical stress is affected by competition. Anyway, the previously mentionned studies show that monospecific bioassays can hardly encompass the true effects of herbicides on natural communities and that environmental risk assessment procedures can thus sometimes not be suitable enough to protect entirely these communities. Those studies however focus on only two species, thus limiting generalisations to other species. This assessment is further complicated by the fact that the applied modelling method usually considers the datasets obtained for the different species and measured endpoints independently. Independent concentration-response curves are built and a critical effect concentration (CEC) as the x percent effect concentration ( EC x ) is derived from each curve to summarize it. Those CEC are finally modelled into species sensitivity distributions (SSD) as for example in [START_REF] Baillard | Effect of interspecific competition on species sensitivity distribution models: Analysis of plant responses to chemical stress[END_REF]. Some works however tried to implement more complex modelling techniques. For example, Kon Kam King et al. implemented a hierarchical model for six herbicides SSDs that include every concentration-response models [START_REF] Kam King | Hierarchical modelling of species sensitivity distribution: development and application to the case of diatoms exposed to several herbicides[END_REF]. One of the main advantages of this approach was that it integrated the whole information of the data in the models and not only a statistical summary as the EC 10 . It thus propagated the uncertainty in the different levels of organization of the data (the different species and contaminant concentrations used) and it gave insights about a global response of the community against the used contaminants. In this work, this approach was carried out further on a dataset previously used in Baillard et al. where herbaceous plants were exposed to isoproturon in presence/absence of a competitor [START_REF] Baillard | Effect of interspecific competition on species sensitivity distribution models: Analysis of plant responses to chemical stress[END_REF]. We indeed tried to assess the effects of competition on plant responses to an herbicide by comparing the entirety of modelled concentration-response curves and analysing if and how the different parameters of said concentration-response curves (the level of the two asymptotes, the steepness of the curve and the position of the inflexion point) were modified by competition. To do this, we modelled at once the concentration-response curves for the different tested species and the two competition modalities to assess in an integrated manner the effects of competition on plant response against chemical stress. Our hypotheses were that, 1) such a modelling framework would be suitable to help understand the interaction between chemical stress and interspecific competition and 2) based on Baillard et al., interspecific competition would negatively affect plant responses to isoproturon [START_REF] Baillard | Effect of interspecific competition on species sensitivity distribution models: Analysis of plant responses to chemical stress[END_REF].

Materials and methods

Experimental design

The dataset used in this work has previously been published in Baillard et al., where further details can be found [START_REF] Baillard | Effect of interspecific competition on species sensitivity distribution models: Analysis of plant responses to chemical stress[END_REF]. Briefly, a series of ecotoxicological bioassays on herbaceous plants were carried out in two situations: 1) with plants isolated in each experimental device (i.e.

without competition) and 2) with plants included in a standardized hexagonal matrix of 37 seedlings of a competitor species Bromus erectus [START_REF] Birch | Rectangular and hexagonal grids used for observation, experiment and simulation in ecology[END_REF]. Six species of the Poaceae family, representative of plant diversity of grassland communities, were used for these experiments: Arrhenatherum elatius, Dactylis glomerata, Lolium multiflorum, Poa pratensis, Poa trivialis and Trisetum flavescens. Microcosms (three-liter-round microcosms, 20 cm diameter, 12.9 cm height containing inert sterilized vermiculite as substrate) were used for the experiments and vegetalized by transferring seedlings at 2-leaves phenological stage, previously grown from seeds in absence of chemical stress. Exposure to isoproturon began after a 4-days acclimatization.

Isoproturon was added to the microcosms via watering to induce chemical exposure to the roots.

The watering protocol was as follows: 150 mL of contaminated nutritive solution [Hoagland basal salt mix (No2, Caisson Laboratories, Smithfield, UT, USA) at 0.82 g/L, pH6] twice a week during the 25 days experiment. Five isoproturon concentrations (0.25, 0.5, 1, 1.5 and 1.75 µM, corresponding to, respectively, 51.5, 103, 206, 309, and 360 µg/L) plus a control were used and for each concentration (six), competition modality (two) and species (six), eight replicates were realized, giving a total of 576 experimental devices.

Endpoint analysis and dataset description

After the exposure, a set of eleven endpoints were measured on each organism. These endpoints were chosen to assess the effects of the two tested factors (interspecific competition and chemical stress) on the species under study. A list of publications used to set up this endpoint list can be found in SI (A). Some endpoints were chosen for their expected response to both competition and herbicide. This was the case of maximal shoot height, root length, root dry mass (DM) and shoot DM, which represent global endpoints of plant performance as an outcome of functional trait responses. Some other endpoints were chosen mainly for their response to competition as functional responsive traits: Leaf Dry Matter Content (LDMC), ligula height (corresponding to stem height) and Specific Leaf Area (SLA). A third group of endpoints was chosen for its response mainly to herbicide: carotenoid and chlorophyll contents and maximum efficiency of photosystem II determined from Fv/Fm (later called Fv/Fm). A final performance endpoint, total DM, was calculated by adding together root and shoot DM. This was considered as a performance trait indicating plant ability to tolerate isoproturon and competition. We also measured B. erectus shoot DM at the end of the experiment to ensure that this DM was constant between the different isoproturon concentrations for each species and that competition pressure was thus constant all along each concentration-response curve.

A log-tranformation was applied on every endpoint but Fv/Fm. For this latter, as it is expressed as a proportion, a logit transformation was applied. Such transformations helped improve data normality and homoscedasticity as shown in [START_REF] Chiarini | Influence of plant development, cultivar and soil type on microbial colonization of maize roots[END_REF] and Baillard et al.

[2020].

Data modelling

Rationale

The main objective was to derive parameter values for the concentration-response curves (CRCs)

and to estimate the modification of those parameter values due to interspecific competition. When looking at CRCs independently, it is impossible to draw a global and statistically significant conclusion concerning the effects of competition. Here, we built a modelling framework allowing to take into account every data simultaneously and to produce, for each endpoint, a global answer to the questions "What is the effect of interspecific competition on herbaceous plant response toward an herbicide?" and "Are the effects of competition similar between the different tested species?". For each endpoint, twelve CRCs were described within one single hierarchical model (six species * two competition modalities), allowing to describe every data point (i.e. for every species, competition modality, isoproturon concentration and replicate) and to study the effects of interspecific competition on the CRCs of all species at once. The modelling framework was built to assess the effects of competition on different elements of the CRC as for example the modification of the level of response at the control, i.e. without chemical stress, or the shift in EC 50 values.

Presentation of the core model

For each endpoint, a modelling framework was designed to encompass twelve CRCs. We used a log-probit model to describe those curves:

      , , log / , , , , , , , , , , , , , , , = * *
= with 0, k i j ij xe i j k i j i j i b i j k l i j k i j k l i j k l i d dc resp yS       (1)
where i indicates the species, j indicates the competition modality (0 means there is no competition and 1 means there is competition), k indicates the isoproturon concentration and l indicates the replicate. given endpoint, so that this parameter is defined at the species level. showing a visible variation of parameters value with competition: the blue curve and the parameters named with the index 0 refers to the situation without competition and the red curve and the parameters named with the index 1 to the situation with interspecific competition.

Modelling of competition effects, parameter definition and model choice

For each species (index i ), two sets of basic log-probit parameters logb (this group of parameters later defined as comp , see SI B). For a given endpoint, we considered that each comp parameter could be of three different types. First, it could be a null effect and thus be equal to 0 for every species as it was possible for interspecific competition to have no visible effects on the parameter for the given endpoint. Second, it could have a fixed value among the different species (generic effect) meaning that modelling the same non-null effect of competition on CRC best described the data. Finally, the competition parameter could have a different value for each species (random effect, i.e. species-specific) as reaction to competition is not necessarily the same for every species. In that last case, competition parameters became respectively , ij  parameters were defined at the species level, meaning that each parameter had a different values for each species, but these values were assumed to follow a same Gaussian law, as frequently hypothesized in Species Sensitivity Distribution modelling [START_REF] Posthuma | Species sensitivity distributions in ecotoxicology[END_REF].

Likewise, when comp parameters were defined as random (i.e. species-specific parameters), they were defined at the species level and assumed to follow Normal laws.

The type of each competition parameter was defined by calculating the DIC (Deviance Information Criterion, a criterion of fit quality for hierarchical models using Bayesian inference penalized by the model complexity) of models built using every different combination types for each competition parameter and selecting the model with the lowest DIC. When DIC values were close (with a difference of less than two, as commonly applied for the AIC criterion in frequentist inferences), the simplest model was kept. For a global view of the built models, a summary of the defined nodes and their definition and an example of a probabilistic Directed Acyclic Graph (DAG) can be found in SI B.

Inference and parameter estimation

Bayesian inference was carried out to estimate model parameters. This method provides estimates as posterior distributions, based on prior distributions and observations. In this work, mainly weakly informative priors based on the experimental design were chosen to keep inferences in a reasonable range. The details of how those priors were defined and the explanations of these values can be found in SI C. Posterior parameter estimation by Bayesian inference were realized with Markov Chains Monte-Carlo techniques [START_REF] Gilks | Introducing markov chain monte carlo[END_REF] using JAGS [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF] and Rjags [START_REF] Plummer | rjags: Bayesian Graphical Models using MCMC[END_REF] under the R environment [R Core [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Three independent chains were generated. For each model, a Raftery and Lewis diagnostic [START_REF] Raftery | How Many Iterations in the Gibbs Sampler?[END_REF] was done to define the number of iterations necessary to reach a good convergence in our simulations. For each endpoint, chain convergence was checked visually and using Gelman and

Rubin convergence diagnostic [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF]. To take into account potential autocorrelation of iterations, we applied a "thin" on the chains, consisting in keeping only one simulated point every thin number. The precise procedure and settings used for simulations can be found in SI D. DIC values were calculated using the same number of iterations and thin. Parameter estimation was realized in two steps since, because of its boolean nature, the estimation of i resp at the same time of the other parameters caused convergence difficulties for these parameters.

First, the simulations were run as described above but only the estimated values of the i resp parameter were kept. Those values were considered to be of good quality after a visual inspection of the CRCs. In a second time, the simulations were run again with i resp values fixed at 0 or 1 according to the results of the first simulation.

Results

Modelling quality indicators

For a majority of endpoints, the procedure was able to return CRCs that described well the data.

The exceptions were for chlorophyll and carotenoid contents that exhibited high variability and for which responses varied very little against the isoproturon concentration. Those endpoints were thus excluded from the analysis (raw data plots for those two endpoints can be found in SI E), thus leaving nine remaining endpoints for the modelling.

With a maximum value of about 1.04, the results obtained for the shows that, for every endpoint, at least one comp parameter was defined to have a fixed or random effect, evidencing that an effect of interspecific competition was visible for each endpoint.

However for two endpoints (root length and root DM) this competition only modified the response at the control ( compd ), meaning that the other aspects of the modelled CECs were not significantly affected by the competition.

Table 1: In the first part: models chosen by the procedure described in this article. It is indicated for each endpoint which was the effect type for each comp or i comp between null (represented by a "-" sign), fixed and random effect. In the second part: summary of the chosen competition effect type of each estimated comp or value of zero, corresponding to a constant curve). Shoot and total DM performance endpoints had a pattern close from the one displayed by SLA. Indeed, competition had an effect on the same two parameters, but here, this effect was random on i compd meaning that competition alone had a different effect on the different species, and a fixed effect on comploge meaning that the horizontal shift of the curves with and without competition was the same from one species to another (SI G). On the four remaining endpoints, three of them, Fv/Fm, LDMC and ligula height, had respectively 1, 2 and 3 parameters with a fixed effect. They had a fixed effect on complogb . In addition, for LDMC and ligula height, a fixed effect was applied for comploge . The comploge and complogb parameters together define the species sensibility against the contaminant and are the most important to assess the effects of competition on species response. Finally, ligula height also had a fixed effect on comp dc  . This parameter is more complicated to interpret as it reflects a change in response amplitude between the two asymptotes of the curves. However, as there was no effect on the compd parameter, it meant that the species were more strongly impacted by high doses of isoproturon when the competition constraint was present. The remaining endpoint is the maximum shoot height on which were applied a random effect for i compd as for shoot and total DM endpoints and a fixed effect on comp dc  . It is difficult to interpret this situation directly. The fitted CRCs for this endpoint (SI G) shows that some species are less affected by high concentrations of isoproturon in presence of a competition constraint. The second part of Table 1 summarizes the number of null, fixed and random effects of competition on each of the four parameters of CRCs for the nine endpoints. We can see in this table that adding a competition effect in the concentration response curves described the data better for 16 parameters on 36. These effects were mainly fixed effects (12 parameters on 36) and random effects were rarely selected by the modelling framework (4 parameters on 36). A small majority of competition parameters was determined as corresponding to null effects (20 parameters on 36), meaning that data were not described better when effects due to interspecific competition was taken into account for these parameters.

Graphs showing the estimated competition effects for each comp and i comp parameter can be found in Figure 3a for fixed effects (i.e. generic effects) and 3b for random effects (i.e. specific effects). We can see on these two figures that observed competition effects could be positive or negative depending on the measured endpoints. In some cases, for example for i compd of shoot DM (Figure 3b), this effect was positive for some species (L. multiflorum) and negative for the others, showing a variety of reaction against competition not only in intensity, but sometimes also in direction. Another important point was that for every random effect on i compd (this trend was less visible on i comploge ), the rank of intensity of this effect between the different species was similar: D. glomerata and P. trivialis showed the most negative effect, A. elatius, P. pratensis and T. flavescens an intermediate effect and L. multiflorum the most positive effect, thus providing information about the difference of response to competitor of these species (species have been ordered on Figure 3b to highlight this point).

Figure 3: Figures showing median (point) and 95% IC (solid line) of the estimated comp and i comp parameter value to the effect type, the endpoint and the species (for random effect only). The vertical grey dashed lines highlight the zero value. For the fixed effect parameters (a), each graphic represents the information for a given comp parameter and each line give the associated value for a specific endpoint. The word "random" means that the parameter was defined to be random and is thus represented on Figure 3b. No information means that the parameter was defined to be of the null type. For the random effect parameters (b), each graphic represents the information for a given i comp parameter and endpoint combination. On each line is shown the informations for a specific species.

Discussion

Despite recent advances [START_REF] Boutin | Effects of Sublethal Doses of Herbicides on the Competitive Interactions Between 2 Nontarget Plants, Centaurea cyanus L. and Silene noctiflora L[END_REF][START_REF] Damgaard | The effect of glyphosate on the growth and competitive effect of perennial grass species in semi-natural grasslands[END_REF], the effects of competition on plant response to chemical stress remain poorly characterised. In this work, we modelled in a single modelling framework data from monospecific bioassays in presence and absence of competitor for six phylogenetically close species. The aim was to see how competition modifies the CRCs of those species for the studied contaminant. This hierarchical modelling showed that an effect of interspecific competition on plant response to herbicide was visible for every species and selected endpoint, highlighting the importance of using multi endpoint and/or multi species approaches in ERA. This effect was applied by the modelling framework on close to a half of the total estimated CRC parameters and only 25% of these effects varied from one species to another, the others being defined as constant between the different species.

This approach had the advantage to estimate the values of every parameters for a given endpoint at once. The CRCs of the six species and two competition modalities are defined all together and are thus interdependent. This is not the case in classical ERA approaches as in Baillard et al. where CRCs were defined independently [START_REF] Baillard | Effect of interspecific competition on species sensitivity distribution models: Analysis of plant responses to chemical stress[END_REF]. In our study, this approach had the advantage to integrate the competition effect for the different species in the four parameters of the sigmoid CRCs, allowing to track down on which aspect the interspecific competition had an effect. An interesting approach would have been to model the whole available dataset inside a unique modelling framework. It was however impossible in those conditions as the different endpoints exhibited a range of values that were each too different from the other, thus limiting the shrinkage possibilities. It would have been possible to transform the data to bring them in the same range, for example by calculating percentage of variation compared to the control, but we made the choice not to superimpose different transformation types as we already applied one to improve residuals normality and homoscedasticity. This would also have been problematic for SLA that varies in the opposed direction against isoproturon compared to the other tested endpoints.

In our modelling process, integrating the competition effect in the models better described the dataset. Indeed, for each endpoint analyzed, at least one parameter varied under competition, which led to a significant effect of competition on half of the 36 parameters tested, each of the nine endpoint being associated to 4 parameters (Table 1). In particular, among the 16 parameters modified by competition, ten corresponded to parameters involved in plant responses to isoproturon ( comp dc  or i comp dc  , comploge or i comploge and complogb or i complogb parameters), whereas six were independent from herbicide exposure ( compd or i compd parameters). In addition to confirm the existence of interspecific competition in the experimental design used, these results thus showed that the presence of competitors modified the way plants respond to chemical stress. Moreover, the responses of shoot-related endpoints to competitors in the absence of chemical stress ( compd or i compd parameters) indicated that competition mainly acted on aerial parts of plants. Indeed, the increase of both shoot height and SLA under competition (Figure 3) is characteristic of response to shading, allowing increasing light interception [START_REF] Fan | Soybean (Glycine max L. Merr.) seedlings response to shading: leaf structure, photosynthesis and proteomic analysis[END_REF], Yang et al., 2019]. This was also relevant with the fact that the experiments were conducted under non-limiting water and mineral nutrient conditions at the root level.

The responses of plant endpoints to chemical stress can depend on the mode of action of the contaminant. Although the existence of multiple targets cannot be ruled out [START_REF] Ramel | Xenobiotic sensing and signalling in higher plants[END_REF], the phenylurea herbicide isoproturon is known to target photosynthesis through PSII, with subsequent generation of triplet chlorophyll, singlet oxygen and ROS cascades [Grouselle et al., 1995, Rutherford and[START_REF] Rutherford | Herbicide-induced oxidative stress in photosystem II[END_REF]. Moreover, root-level application of herbicides can lead to shoot impact by whole-plant distribution through the xylem flow [START_REF] Sulmon | Sucrose amendment enhances phytoaccumulation of the herbicide atrazine in Arabidopsis thaliana[END_REF].

Given its effects on PSII, the present isoproturon exposure can therefore be seen as a situation of leaf PSII perturbation and photooxidative stress with a cascade of regulatory and developmental consequences involving metabolic, hormonal and shoot-root signalling mechanisms. Aerial parts of the plants were thus directly targeted by both isoproturon treatment and interspecific competition (through shading). This explains the results of Table 1 where isoproturon-related parameters exhibited a competition effect only for the shoot-related endpoints. In the model, the , ij e parameter corresponds to the EC 50 value of the concentration-response curve. The comploge or i comploge parameters, estimating the modification of CRC EC 50 values under competition, can thus be taken as a relevant marker of chemical stresscompetition interaction. These parameters were found to be significant for the shoot-related endpoints LDMC, SLA, ligula height, and shoot and total DM (Table 1), competition inducing in each case a decrease of EC 50 values (Figure 3, SI G). As functional traits directly involved in competition response [START_REF] Fan | Soybean (Glycine max L. Merr.) seedlings response to shading: leaf structure, photosynthesis and proteomic analysis[END_REF], Yang et al., 2019], the leaf-dedicated LDMC and SLA endpoints were particularly affected by competition. This makes them particularly relevant endpoints to conclude that interspecific competition increases the inhibition effect of isoproturon on plant species or makes them more sensitive to isoproturon. In contrast, shoot and total DM, which also exhibited a competition effect on , ij loge , consisted in more global endpoints of plant performance, resulting from the integration of functional trait responses to competition and chemical stress exposure. The complogb and i complogb parameters also demonstrated interesting effects of interspecific competition on mainly leaf-related endpoints (Table 1). Competition was thus found to increase CRC steepness of Fv/Fm whereas decreasing it for LDMC (SI G). These results thus suggest that competition relieves isoproturon effect on photosynthesis, as herbicide-targeted process, and emphasizes its effect on a more global leaf physiology-related endpoint (LDMC). Competition, by inducing leaf shading, has been shown to change leaf physiology, leading to an increase of chloroplast number and grana layer thickness [START_REF] Fan | Soybean (Glycine max L. Merr.) seedlings response to shading: leaf structure, photosynthesis and proteomic analysis[END_REF]. As these changes may increase the number of PSII and consequently the number of D1 proteins targeted by isoproturon, the herbicide would thus be less efficient in inhibiting a higher number of D1 targets. LDMC represents a proxy of leaf density [START_REF] Ryser | The Importance of Tissue Density for Growth and Life Span of Leaves and Roots: A Comparison of Five Ecologically Contrasting Grasses[END_REF], which is strongly related to its biochemical composition and in particular to its protein content [START_REF] Wilson | Specific leaf area and leaf dry matter content as alternative predictors of plant strategies[END_REF]]. Yet, this endpoint has been shown to be decreased by shading under competition [START_REF] Yang | Effects of shading on the growth and leaf photosynthetic characteristics of three forages in an apple orchard on the Loess Plateau of eastern Gansu, China[END_REF]. Thus, changes in leaf components, and especially the decrease in leaf protein content suggested by LDMC decrease under competition, may explain the increase of plant sensitivity to isoproturon since plants need numerous cell compounds (antioxidants, chaperone proteins, antioxidant enzymes) to respond to chemical stress and particularly to PSII inhibiting herbicides [START_REF] Ramel | Differential patterns of reactive oxygen species and antioxidative mechanisms during atrazine injury and sucrose-induced tolerance in Arabidopsis thaliana plantlets[END_REF][START_REF] Ramel | Genome-wide interacting effects of sucrose and herbicide-mediated stress in Arabidopsis thaliana: novel insights into atrazine toxicity and sucrose-induced tolerance[END_REF]. The , ij dc  parameter was the least affected by competition. Indeed, this parameter, which estimates the inhibitory effect of the highest isoproturon treatment (1.75µM), exhibited a positive competition effect only for the shoot-related height endpoints (ligula height, max shoot height; Table 1, Figure 3a), which were more inhibited by isoproturon under competition than without. It thus seemed that competition weakly affected plant responses at high chemical stress intensity, in accordance with previous studies on abiotic stresses demonstrating that competition acts mainly under low stress intensity [START_REF] Callaway | Competition and Facilitation: A Synthetic Approach to Interactions in Plant Communities[END_REF][START_REF] Choler | Competition and facilitation on gradients in alpine communities[END_REF][START_REF] Jung | Impact of competition on plant performances along a flooding gradient: a multi-species experiment[END_REF][START_REF] Liancourt | Stress Tolerance and Competitive-Response Ability Determine the Outcome of Biotic Interactions[END_REF].

Moreover, the fact that competition had reversed effects on isoproturon response for shoot height, adverse for low concentrations and advantageous for high concentrations (SI G) strongly suggested the involvement of energetic trade-offs under combined chemical and biotic stresses.

This study was carried out by using six species of the Poaceae family in order to be representative of some plant biodiversity in grass communities and in order to consider the potential diversity of responses to chemical stress and to competitors. This experimental design thus showed that most competition effects identified for parameters related to chemical stress responses (i.e. effects on , ij dc  ,

, ij loge and , ij logb ) were characterized as fixed (i.e. as comp dc  , comploge and complogb ), and thus species-independent. In contrast, three out of the four random competition effects were applied on , ij d parameter, which is chemical stress-unrelated. For most of the endpoints analyzed, it thus seemed that, under chemical stress, the competition effect remained constant whatever the species targeted, and whatever its effect in the absence of chemical stress, thus suggesting a role of chemical stress in homogenizing competition effects. Given the focus of our study on Poaceae, the characterization of generic parameters of chemical stress-competition interactions will be particularly useful for ERA and plant dynamics modelling of grazed pasture or grassland plant communities, which cover a large fraction of the Earth's land surface [START_REF] Barbehenn | C3 grasses have higher nutritional quality than C4 grasses under ambient and elevated atmospheric CO2[END_REF]. Moreover, induction of situation of PSII perturbation is also known to occur in a number of abiotic stresses, such as high light, drought, salinity, or heat [START_REF] Pospíšil | Formation of singlet oxygen and protection against its oxidative damage in Photosystem II under abiotic stress[END_REF], which are known to be significant components of climate change constraints [START_REF] Bigot | Pivotal roles of environmental sensing and signaling mechanisms in plant responses to climate change[END_REF]. Competition-herbicide interactions at the level of PSII therefore also suggest the potential interference of competition with climate change responses [START_REF] Bigot | Pivotal roles of environmental sensing and signaling mechanisms in plant responses to climate change[END_REF][START_REF] D'alessandro | Sensing  -carotene oxidation in photosystem II to master plant stress tolerance[END_REF][START_REF] Munné-Bosch | The Impact of Global Change Factors on Redox Signaling Underpinning Stress Tolerance[END_REF][START_REF] Noctor | Climate Change, CO2, and Defense: The Metabolic, Redox, and Signaling Perspectives[END_REF].

Conclusion

This study helped understand better how competition modifies the way that herbaceous plants respond to chemical stress and how this response varies from one species to another. Competition effects were observable on every endpoint and these effects could be tracked down on the different parameters of the fitted concentration-response curves. This innovative approach had the major advantage to integrate every information given by the dataset and to propagate uncertainty from one hierarchical level to another. This work thus gives new insights for the integration of realistic supplementary factors, such as competition, into ecotoxicological bioassays and environmental risk assessments.

Graphical abstract

Highlights

In ecotoxicology, dose-response curves are used as species-and condition-specific We modelled curves for six grassland plant species and two competition modalities Hierarchical modelling integrated competition to species and contaminant effects Competition effect was significant and mostly generic rather than species-specific Modelling is a way to handle multiple factor effects in ecotoxicological bioassays

  standard deviation of the law of concentrations in log scale of which  is the cumulative distribution function (assumed Gaussian). Hence, the lower is , ij b , the steeper is the log-probit model. A summary of the role of these different parameters on the log-probit curves can be found in Figure 1used hereafter. For each species, replicate spreading around CRC is considered to be Gaussian and is characterized by a standard deviation i S supposed to be constant inside a given species and thus independent of concentration and competition. The i resp parameter indicates whether the log-probit model is necessary to describe the concentration-response relationship or if a constant model at the , ij d level describes the data sufficiently well. We considered that i resp should be the same with and without competition for a

Figure 1 :

 1 Figure 1: Scheme of the role of the different parameters of the basic log-probit model, their

  group of parameters later defined as i comp parameters) as they had different values for the different species. The effects of those competition parameters on the CRC can be highlighted and translated into biological terms. compd translates the modification of the response level at the control (i.e. without chemical stress) and thus shows the effects of competition alone on the monitored endpoint. The three other comp parameters ( comp dc  , comploge and complogb ) are for their part translating a modification of the toxic effects by competition ( being contaminant-related parameters) as comp dc  translates a change in the amplitude of the toxic response, comploge a change in the concentration at which is situated the inflection point of the curve (and thus the effects of competition on EC 50 which is commonly used in ERA) and complogb a change in the concentration-response curve steepness.

  Figure2aand for SLA in Figure2b(the CRCs fitted for every other endpoint are provided in SI G).These two examples were chosen because they were good examples of how the different effect types applied by the modelling framework were translated in CRCs. For root DM (and for root length), the data were better described with a null effect of competition on

Figure 2 :

 2 Figure 2: Example of concentration-response curves produced using parameters of the model built for the Root DM and SLA endpoints. Open circles represent the data point, plain circles represent the mean of data points for each concentration, the lines are the CRCs drawn using the mean estimated value of each parameter and the vertical plain lines are the 95% Credible Interval (CI) of the model response at each tested concentration (abscissa values were modified for those lines not to superimpose on the data points). In blue are the data/model without competition and in red the data/model with competition.

Figure 1 Figure 2 Figure 3

 123 Figure1

  Gelman and Rubin indicator were low enough to indicate a good convergence of the Bayesian inference. Detailed results are provided in SI F. Table 1 shows which effect type was applied on each ,0

	i  parameter
	(in absence of competition) to obtain its ,1

i  counterpart (in presence of competitors). This table

Acknowledgements

This work has been funded by the "Ecosphère continentale et côtière" (EC2CO) interdisciplinary program from the Centre National de la Recherche Scientifique (CNRS, France). Financial support for mobility has also been provided by the Groupe de Recherche (GDR) 3645 Statistical Ecology (CNRS, France) to Vincent BAILLARD. Vincent BAILLARD finally received financial support for salary from French Research ministry and the Région Grand Est. The authors are grateful to Valerie GOUESBET (Ecobio Research Unit, Experimental Ecology platform), to Lucie WRONKA and Malo BOURGET for help with experimentations, and to Sandra RIGAUD (Ecobio Research Unit) for administrative and accounting work. The authors thank Jean-Pierre CAUDAL (Ecobio Research Unit, Experimental Ecology platform) and the Experimental Ecology platform of Ecobio research unit for providing all the infrastructure needed to carry out the experiments.

Declaration of competing interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: