Joel Jose Montalvo Acosta

Dr Marian Dryzhakov

Dr Edward Richmond

Dr Marco Cecchini
email: mcecchini@unistra.fr

Prof. Dr Joseph Moran
email: moran@unistra.fr

Dr Joel J Montalvo-Acosta

A Supramolecular Model for the Co-Catalytic Role of Nitro Compounds in Brønsted Acid Catalyzed Reactions

Keywords: Brønsted acid catalysis, molecular modeling, nitro compounds, H-bond aggregation, rational de novo design

come

The proton provides the most common and important way to catalyze organic reactions, [1] however it needs not always operate in isolation. Within a single molecular catalyst, such as in the active site of an enzyme, multiple hydrogen-bond donors (or Brønsted acidic sites) can work in cooperation to increase the overall acidity and hence catalytic activity. [2] Cooperative Hbonding and Brønsted acidity is also possible between molecules, when multiple H-bond donors or Brønsted acid molecules interact to generate an aggregate that is a more effective catalyst than either individual molecule. [3] However, what role the bystander molecules -those that are not H-bond donors or Brønsted acids -might play, if any, remains to be elucidated. Though poorly understood, a few papers indicate that interactions between Brønsted acids and seemingly innocuous solvents or additives, particularly nitro compounds, can have a profound influence on reactivity both in terms of accelerating reaction rates and by changing the concentration dependence of the reaction. In the 1960s, Pocker and coworkers found that the hydrochlorination of olefins was not only markedly faster when carried out in nitromethane compared to other solvents, but also displayed an atypical second order kinetic concentration dependence on HCl. [4] Further investigations to uncover why this effect would be exclusive to nitromethane were not made. [START_REF] Sergeev | Some have attributed the unusual kinetic order to the polarity of the bulk solvent, though no explanation is given as to why other polar solvents do not induce a similar effect[END_REF] In 2015, some of us observed that catalytic quantities of nitro compounds in benzene were highly specific to accelerate the dehydroazidation of tertiary aliphatic alcohols catalyzed by B(C6F5)3•H2O (BCF), [6] a strong Brønsted acid of comparable strength to HCl. [START_REF] Danopoulos | C6F5)3•H2O has a pKa of 8.4 in MeCN and has been studied in detail[END_REF] The reaction was found to display a second order concentration dependence with respect to the nitro compound as well as a second order concentration dependence with respect to BCF. In contrast, when carried out in the absence of nitro compounds, the reaction was much slower and found to be first order with respect to BCF (Figure 1). Snyder ' s group has also observed Brønsted acid-catalyzed reactivity that is highly specific to wet nitromethane. [8] No model currently exists to account for why the presence of a nitro compound, sometimes present only in catalytic quantities, could change the concentration dependence and rate of at least these two different Brønsted acid catalyzed reactions. Here, using DFT calculations and vibrational analysis, we obtain a structural model for a higher order aggregate formed by nitro compounds and Brønsted acids, examine the effect of aggregation on acidity, show that the model can be used to predict reaction rates for a set of new nitro compounds in alcohol dehydroazidation, and even predict a new class of "template" molecules that can induce a similar co-catalytic effect. Beyond the specific reaction in question, the present results imply that a comprehensive understanding of Brønsted acid catalyzed reactions can only be achieved by taking into account the weak interactions between Brønsted acids and all molecules in solution. The dehydroazidation of tertiary alcohols exhibits a secondorder kinetic dependence on the concentration of BCF and a second-order kinetic dependence on nitro compound. [6] As there exists significant spectroscopic and catalytic evidence that nitro compounds can act as H-bond acceptors, [9] we postulated that the reaction kinetics could be explained by the intermediacy of a catalytic competent aggregate involving H-bonded self-assembly of two BCF with two nitro-compound molecules. Structural models of the tetrameric (2:2) assembly of BCF with nitromethane were generated by connecting the molecules via pairs of H-bonds and optimizing the geometry of the aggregate by DFT. Given the complexity of the configurational space to be explored, an initial conformational search was carried out using the HF-3c semi-empirical method, [START_REF] Sure | Corrected small basis set Hartree-Fock method for large systems[END_REF] then refined by DFT calculations using the ωB97X-D functional [11] with the 6-31G(d,p) basis set; see SI for details. DFT optimization results in an almost flat, rectangular hydrogen-bonded network formed between the nitro groups of two nitromethanes and two boron hydrates (Figure 2). The tetrameric arrangement is stabilized by two sets of non-equivalent H-bonds, two shorter and two longer.

A Hydrochlorination of olefins (ref. 4) R 4 R 3 R 1 R 2 HCl solvent H Cl R 1 R 2 R 3 R 4 B Dehydroazidation of tertiary alcohols (ref. 6) OH R 3 R 1 R 2 B(C 6 F 5 ) 3 •H 2 O (cat.) TMSN 3 (3 equiv) N 3 R 3 R 1 R 2 C 6 H 6 ,
Vibrational analysis of the DFT-optimized supramolecular structure indicates that each pair of O-H groups in the network present two distinct stretching modes (i.e. one symmetric with a lower IR intensity signal and one anti-symmetric with higher intensity) with the symmetric stretching of the hydrogens involved in the shorter H-bonds characterized by the lowest vibrational frequency of 3501 cm -1 ; see Figure S1. The same analysis of one isolated BCF or its 1:1 complex with nitromethane shows that the corresponding O-H stretching modes have vibrational frequencies of 3693 and 3670 cm -1 , respectively (Figure S2-S3). Therefore, the formation of a tetrameric adduct with nitromethane results in a red-shift of the O-H stretching frequency as large as 192 cm -1 relative to isolated BCF and 169 cm -1 relative to the 1:1 complex. Since an O-H stretching frequency reduction was shown to provide a measure of the strength of a Brønsted acid in substituted p-and m-phenols (i.e. the lower the frequency, the stronger the acid), [12] our calculations suggest that the 2:2 BCF:nitromethane selfassembly increases the acidity of BCF, in particular for the hydrogens involved in the shorter hydrogen bonds. This counterintuitive result implies that aggregation of multiple acid molecules templated by very weak H-bond acceptors may have the net effect of increasing the acidity, and hence catalytic activity, of the aggregate relative to the single acid molecules. Following the same modeling procedure, 3D structures of the tetrameric assembly of BCF with a variety of nitro compounds including 4-nitrobenzotrifluoride, nitrobenzene, 4-nitroanisole, 2methyl-2-nitropropane and 1-nitrohexane were generated; these compounds were previously shown to modulate the co-catalytic activity of BCF. [6] Upon geometry optimization at the DFT level of theory, the rectangular hydrogen-bonding network made of pairs of non-equivalent H-bonds is largely preserved in all complexes, with minor distortions in some cases; see Figure S5. Vibrational analysis of the model complexes shows the same characteristic red shift in the O-H stretching frequency with a magnitude that is dependent on the chemical nature of the nitro compound; see Table S1. Remarkably, the red shift in the anti-symmetric stretching of the short H-bonds is found to be strongly correlated with the experimental Log(rate) with a determination coefficient (R 2 ) of 0.81; see Figure 3 (filled circles). This observation suggests that the O-H stretching frequencies predicted by DFT in the 2:2 aggregate in vacuum can be used as predictors for the kinetic rate of dehydroazidation catalyzed by BCF. The same calculations in benzene for a subset of compounds using an implicit-solvent model indicate that solvation effects in apolar media play only a minor role on the structure and the O-H stretching frequency of the 2:2 assembly thus justifying the use of gas phase calculations; see SI. In addition, an analysis of the level of theory on the quantification of the O-H stretching frequency shift upon self-assembly indicates that the use of the 6-31G (d,p) basis set provides reliable predictions at the lowest computational effort; see SI.

To evaluate the statistical significance of the correlation in Figure 3, two validation schemes were followed: [13] a cross internal validation based on a leave-one-out (LOO) analysis; and an external validation based on a new set of five nitro compounds (1-nitropropane, 2-nitropropane, nitroethane, 2nitroanisole and 1-nitronaphtalene), whose influence on reaction rate was previously unknown and determined by GC-MS; see SI. The internal validation indicates that the DFT model presents a positive predictive profile with a cross-validated squared correlation coefficient (q 2 ) of 0.46 (Table S3). The external validation shows that kinetic rates predicted by the DFT model are strongly correlated with the experimental rates with R 2 of 0.70; see Figure 3 (empty circles). In addition, the statistical parameters of the DFT model appear to be robust to randomization in alternative training/test sets, which is consistent with the absence of random correlations (Table S4). The predictive character of the DFT model for the experimental Log(rate) supports the assertion that tetrameric self-assembly of BCF with nitro compounds is a critical feature of the catalytic mechanism in the azidation of tertiary alcohols. [6] A similar model for the Log(rate) was constructed using experimentally observed IR stretches, i.e. the kinetic rate of azidation was correlated with the frequency at the maximum of the OH-stretching band in the experimental IR spectrum; see SI for details. Although this IR model was predictive for the initial set of nitro compounds, it was not for the test set (Figure S7). In particular, the accelerations produced by 2-nitroanisole and 1nitronapthalene (entries 11-12 in Table S1) are reversed by the IR model, which produces an anti-correlation between predicted hb4 hb1 hb2 hb3 and observed Log(rate) in the external validation; see Figure S8. The more limited predictive power of the IR model compared to DFT is explained on the basis that the intense OH-stretching bands derive from the overlap of multiple stretching modes, which limits the resolution of the experimental determination. Despite the inaccuracies associated with the DFT functional/basis set [14] , the computational vibrational analysis uniquely allows to quantify the red shift in the anti-symmetric stretching of the short H-bonds, which is predictive for the Log(rate); see Figure 3. S1.

The formation of the catalytically active species isolated by calculations involves a largely unfavorable entropic penalty due to the association of four molecular entities. To test the thermodynamic stability of the 2:2 assemblies, the free energy of association was investigated in the rigid-rotor harmonic oscillator (RRHO) approximation; see SI. The results show that all 2:2 complexes involving nitro compounds are marginally stable with standard association free energies ranging from +3 to +14 kcal/mol; see Table S7. Although the 2:2 complexes are predicted as rarely populated in solution, we reason that if the self-assembly were faster than the background reaction, one would still observe a co-catalytic effect. Based on this hypothesis and using the free energy of self-assembly as a lower-bound estimate for the kinetic barrier of tetrameric association, we evaluated the activation free energy of the background reaction using Eyring's equation [START_REF]The activation free energy for the background reaction, i.e. dehydroazidation in the absence of nitro compounds, was estimated from the experimental Log(rate) via Eyring's equation k = (kBT/h)•exp(-ΔG° ‡/RT). Using Log(k) = -5.2 (compound 7 in Table S1), Eyring's formulation predicts an activation free energy at 298[END_REF] . The results show that in the absence of nitro compounds the activation barrier for dehydroazidation is at least 10 kcal/mol higher than the kinetic barrier for self-assembly. We suggest that albeit rarely populated, 2:2 self-assembly may still have a detectable effect on the dehydroazidation rate. Last, the DFT model was used to explore new chemotypes capable of playing the same co-catalytic role as nitro compounds. For this purpose, six representative compounds were selected on the basis of their ability to weakly accept hydrogen bonds on two different atoms (Table 1). For each compound, the tetrameric assembly with BCF was modeled and the relevant O-H stretching frequencies calculated by DFT. In four out of six cases, the assembly preserved a rectangularshape and a planar hydrogen-bonded network after geometry optimization. Of these four, the Log(rate) predicted by the DFT model was in good correspondence with the observed reactivity after 24 h under standard reaction conditions in three cases (entries 1-3). Also, the two compounds that did not result in stable tetrameric complexes in the calculations did not show any accelerating effect on the reaction (entries 5-6). Most importantly, as predicted by the DFT model, diethylsulfate (entry 1) as additive showed a similar accelerating effect to nitromethane; [START_REF]activating agent or that trace acidic impurities were responsible for catalysis[END_REF] the tetrameric assembly computed for diethylsulfate is shown in Figure 4. Indeed, the method of initial rates showed a kinetic order dependence of 2.55 on diethylsulfate, consistent with higher order aggregates being involved in catalysis; see SI. ---<5%

[a] A complex was considered stable if it preserves a rectangular-shape and planar hydrogen bond network as obtained in tetramer self-assembly of BCF with nitro compounds. [b] Difference between the OH-stretching frequency of BCF in the complex and alone (3692.52 cm -1 ). [c] Reaction monitored by GCMS relative to dodecane as internal standard.

Notably, this is the first time that an additive other than a nitro compound has been shown to induce changes in the kinetic order dependence and an accelerating effect in the azidation reaction. The incorrect prediction for Log(rate) in entry 4 illustrates that above a certain threshold of H-bond accepting ability, the buffering effect of the additive dominates any positive effects arising from aggregation; see SI for a discussion. In summary, DFT modeling of 2:2 H-bonded aggregates of nitro compounds and BCF yields O-H stretching frequencies that are predictive of the experimentally observed rates in reactions that use the corresponding nitro compound as promoter and can be used to guide the rational design of new co-catalysts. The structural model of the aggregate was used to identify sulfate esters as a new class of promoters, which was verified experimentally a posteriori. The correlation between experiments and calculations suggests that the increased acidity of such aggregates could indeed be responsible for the higher order kinetic concentration dependence observed in Brønsted acid catalyzed alcohol dehydroazidation. [6] It remains to be determined whether this supramolecular mechanistic model can also be extrapolated to other reactions shown to have higher order dependence on Brønsted acids when nitro compounds are present, such as olefin hydrochlorination. [4] More broadly, the important implication of this model is that a deep understanding of Brønsted acid catalysis requires consideration of not only the molecular structure and pKa of the Brønsted acid or the bulk properties of the solvent, but also the microscopic nature of the supramolecular environment that takes into account the weak interactions between all molecules in the reaction mixture.

Entry for the Table of Contents

Catalysis by aggregates could explain the unusual higher order kinetics in certain Brønsted acid catalyzed reactions. Computed IR stretches of the aggregate predict reaction rates and new classes of co-catalysts (sulfate esters).

! Catalysis by aggregates could explain the unusual higher order kinetics in certain Brønsted acid catalyzed reactions ! Computed IR stretch of the aggregate predicts observed reaction rates and new co-catalyst structures (sulfate esters) sulfate esters nitro compounds

Figure 1 .

 1 Figure 1. The effect of nitromethane on the kinetic concentration dependence of reactions consuming (A) or catalyzed (B) by Brønsted acids.

Figure 2 .

 2 Figure 2. DFT-optimized structure for the 2:2 tetrameric self-assembly of BCF and nitromethane. The calculated structure highlights the existence of two sets of non-equivalent H-bonds, two shorter ones (hb1 and hb3) and two longer ones (hb2 and hb4).

Figure 3 .

 3 Figure3. Correlation of the experimental kinetic rate with the DFT calculated frequency of the anti-symmetric stretching of the short H-bonds in the 2:2 aggregate. Filled and empty data points correspond to the training (six + blank) and test (five) sets of nitro compounds investigated here. The numbering is the same as in TableS1.

Figure 4 .

 4 Figure 4. Optimized structure for the complex BCF:diethylsulfate (2:2).

Table 1 .

 1 Calculated frequency analysis for tetrameric self-assembly of BCF with various promoters and comparison to the observed reactivity in the dehydroazidation of alcohols.

	entry promoter	OH stretch	b Δν	predicted	conv. at 24
	(stable	(cm -1 )	(cm -1 )	Log(rate)	h (%) c
	aggregate? a )
	1	3510.47	182.05	-4.66	>99%
	(yes)
	2	3561.43	131.09	-4.85	60%
	(yes)
	3	3562.81	129.71	-4.85	58%
	(yes)
	4	3449.12	243.4	-4.44	<5%
	(yes)

Acknowledgements

This work was supported in part by the Agence National de la Recherche through the LabEx "Chemistry of Complex Systems" (ANR-10-LABX-0026 CSC; CSC-MCE-15). Financial support from the International Center for Frontier Research in Chemistry (icFRC) is gratefully acknowledged.


Florent Dewez

Benjamin Guedj

Arthur Talpaert

Vincent Vandewalle

An end-to-end data-driven optimisation framework for constrained trajectories

Keywords: Statistical modelling, Functional data, Constrained optimisation

Many real-world problems require to optimise trajectories under constraints. Classical approaches are based on optimal control methods but require an exact knowledge of the underlying dynamics, which could be challenging or even out of reach. In this paper, we leverage data-driven approaches to design a new end-to-end framework which is dynamics-free for optimised and realistic trajectories. We first decompose the trajectories on function basis, trading the initial infinite dimension problem on a multivariate functional space for a parameter optimisation problem. A maximum a posteriori approach which incorporates information from data is used to obtain a new optimisation problem which is regularised. The penalised term focuses the search on a region centered on data and includes estimated linear constraints in the problem. We apply our data-driven approach to two settings in aeronautics and sailing routes optimisation, yielding commanding results. The developed approach has been implemented in the Python library PyRotor.

Introduction

Many real-world problems require to optimise trajectories under constraints. The present paper stems from an initial work on aeronautics, and the quest for designing fuel efficient aircraft trajectories based on available flight data. We have reached a generic data-driven methodology which falls in the much broader field of trajectory optimisation under constraints. As such, it has potential applications to many real world problems, such as in robotics to minimise the work-based specific mechanical cost of transport [START_REF] Srinivasan | Computer optimization of a minimal biped model discovers walking and running[END_REF] or in aerospace to reduce the total thermal flux when a space shuttle re-enters in the atmosphere [START_REF] Trélat | Optimal control and applications to aerospace: Some results and challenges[END_REF].

In aeronautics, optimisation problems are often formulated in terms of optimal control problems [START_REF] Dalmau | How much fuel and time can be saved in a perfect flight trajectory ? continuous cruise climbs vs. conventional operations[END_REF][START_REF] Girardet | Wind-optimal path planning: Application to aircraft trajectories[END_REF][START_REF] Cots | Direct and indirect methods in optimal control with state constraints and the climbing trajectory of an aircraft[END_REF]. They can be solved by converting the problem into a parameter optimisation problem. This allows to take into account the dynamics of the system, leading to realistic solutions complying with additional constraints (we refer to Rao, 2009, for an overview).

Nevertheless the differential equations describing the dynamics of the system of interest may be (partially) unknown. For instance, the differential system describing the motion of an aircraft moving in an air mass [START_REF] Rommel | Gaussian mixture penalty for trajectory optimization problems[END_REF] involves the lift and drag forces for which no analytic formulas exist. Aircraft manufacturer computes numerical models by means of heavy simulations and wind tunnel tests. Another approach consists in reconstructing unknown forces based on physical formulas and available flight data; see for instance [START_REF] Rommel | Aircraft dynamics identification for optimal control[END_REF]; [START_REF] Dewez | From industry-wide parameters to aircraft-centric on-flight inference: improving aeronautics performance prediction with machine learning[END_REF] for results in aeronautics and [START_REF] Ramsay | Parameter estimation for differential equations: a generalized smoothing approach[END_REF] in the generic setting of parameter estimation for differential equations. However, while being promising, this reconstruction step requires restrictive assumptions and the statistical errors may impact strongly the solution of the optimal control problem. Moreover it does not tackle directly the optimisation problem.

The above approaches require intensive computations and may be affected by noise. The aim of our work is to provide realistic trajectories without involving complex and noisy dynamical systems. It is flexible enough and easily interpretable by experts while permitting the use of efficient optimisation algorithms. In the present work, we leverage available trajectory data to propose a thorough methodology which fulfils the above requirements. Our approach learns a model based on the observed trajectories, allowing in particular to extract some intrinsic properties in a data-driven way. It incorporates this modelling into an optimisation problem through a Bayesian approach. The resulting problem turns out to be constrained by the data in a simple and natural way. The main benefit on this approach is that it directly uses the information contained in the data, requiring no explicit information on the dynamics.

The methodology presented in this paper is specific to the situation where the user has access to trajectory data but, at the same time, the approach is intended to be generic enough so that it can be exploited in a wide range of applications. In particular it is certainly not restricted to the aeronautic setting.

Our approach first assumes that all the trajectories belong to a finite-dimensional space, which allows to reduce the complexity of the problem with low information loss for a wellchosen basis. In a Bayesian framework, we assume that the prior distribution of trajectories (through their related coefficients) is proportional to a decreasing exponential function of the cost, assuming that efficient trajectories are a priori more likely than inefficient ones. The observed trajectories, that we call reference trajectories, are interpreted as noisy observations of an efficient one, the noise following a centered Gaussian multivariate distribution. In a Bayesian perspective it is thus possible to deduce the posterior distribution of the efficient trajectory given the reference trajectories. For the sake of simplicity, we focus on the mode of the posterior distribution, the related objective function involves the cost of a trajectory and its squared Mahalanobis distance to a weighted average of reference trajectories. It can be interpreted as a penalised optimisation problem.

The role of the penalisation in the objective function is to force the solution to be close to real trajectories. The strength of the penalisation is here controlled by a hyper-parameter and a tuning process is proposed to find an optimal balance between optimisation and (non-linear) constraints verification. Hence the optimised trajectory may inherit a realistic behaviour from the above closeness, even though the dynamics are not taken into account in our problem. Further we note that this penalised term is actually quadratic and the optimisation problem is constrained by affine functions. So in certain cases, the problem is convex allowing to make use of very efficient algorithms.

A last remark on the penalised term is the fact that the underlying metrics is given by a covariance matrix estimated on the reference trajectories. It is noteworthy that this matrix not only indicates the most unconstrained directions for the optimisation but also reveals linear relations between variables. In particular, some of these relations may reflect the dynamics or may not be known by the user.

In a nutshell, this data-driven approach restricts the search space to a region centred on the data in a metric space reflecting features estimated from the data. In particular this property may help non-linear optimisation solvers to converge with a limited number of iterations.

Outline. We first describe our approach in Sec. 2. A brief description of the Python library we have develop is provided in Sec. 3. Sec. 4 and Sec. 5 are devoted to applications: the first one to the fuel reduction of aircraft during the climb and the second one to the maximisation of the work of a force field along a path. We finish the paper by discussing on future works to improve and generalise our optimisation methodology.

An end-to-end optimisation workflow based on observed trajectories

We are interested in finding a trajectory y which minimises a certain cost function F , namely a solution of the following optimisation problem:

y ∈ arg min y∈A G (y 0 ,y T ) F (y) . (1)
The set A G (y 0 , y T ), which is defined in Sec. 2.1, models the constraints the trajectory has to comply with, such that the initial and final conditions. Note that a trajectory is typically a multivariate function defined on an interval and its components are given by states and controls (which are not distinguished in this paper for the sake of presentation). If the constraints do not include the dynamics then the solution y may be by far unrealistic. A strategy to force a more realistic pattern would be to add more user-defined constraints to the problem (1), although this may be a complicated task and solving numerically the resulting problem could be computationally expensive.

In view of this, we provide in this section our full workflow to obtain a new optimisation problem which includes in a natural and simple way constraints coming from the data. This problem is actually designed to provide trajectories which have a realistic behaviour.

We begin with elementary but necessary definitions for trajectories and constraints in Sec. 2.1. We aim at stating the optimisation problem in a finite basis space so we define in Sec. 2.2 the mathematical formalisation of how we decompose each trajectory as a projection on such a space. To extract information from the data for the optimisation problem, a statistical modelling on the projected space of the available trajectory data is done in Sec. 2.3. In Sec. 2.4, we put everything together to obtain our new optimisation problem Sec. 2.4 via a maximum a posteriori approach. Sec. 2.5 presents a handy computation regarding the cost function in a quadratic case, for the sake of completeness. Additional details can be found in the supplementary material. Sec. 2.6 focuses on a hyperparameter tuning for an optimal tradeoff between optimisation and additional (nonlinear) constraints. Last but not least, Sec. 2.7 contains confidence intervals to assess the accuracy of the predicted optimised cost when the cost function is known up to a random noise term.

Admissible trajectories modelling

We start with definitions.

Definition 2.1 (Trajectory). Let T > 0 be a real number and let D 1 be an integer. Any continuous R D -valued map y defined on [0, T ], i.e. y ∈ C [0, T ], R D , is called a trajectory over the time interval [0, T ]. The d-th component of a trajectory y will be denoted by y (d) . As such, a trajectory is at least a continuous map on a finite interval.

When optimising a trajectory with respect to a given criterion, the initial and final states are often constrained, that is to say the optimisation is performed in an affine subspace modelling these endpoints conditions. This subspace is introduced just below. Definition 2.2 (Endpoints conditions). Let y 0 , y T ∈ R D . We define the set D(y 0 , y T ) ⊂ C [0, T ], R D as follows:

y ∈ D(y 0 , y T ) ⇐⇒ y(0) = y 0 y(T ) = y T .
In many applications, the trajectories have to satisfy some additional constraints defined by a set of (nonlinear) functions. For instance these functions may model physical or userdefined constraints. In this paper, this set is not intended to include the dynamics of the system. We define now the set of trajectories verifying such additional constraints.

Definition 2.3 (Additional constraints). For = 1, . . . , L, let g be a real-valued function defined on R D . We define the set G ⊂ C [0, T ], R D as the set of trajectories over [0, T ] satisfying the following L inequality constraints given by the functions g , i.e.

y ∈ G ⇐⇒ ∀ = 1, . . . , L ∀ t ∈ [0, T ] g y(t) 0 .
To finish we introduce the set of admissible trajectories which satisfy both the endpoints conditions and the additional constraints.

Definition 2.4 (Admissible trajectory). We define the set A G (y 0 , y T ) ⊂ C [0, T ], R D as follows:

A G (y 0 , y T ) := D(y 0 , y T ) ∩ G .

Any element of A G (y 0 , y T ) will be called an admissible trajectory.

Projection for a finite-dimensional optimisation problem

In our approach, a theoretical optimisation problem in a finite-dimensional space is desired to reduce the inherent complexity of the problem. This can be achieved by decomposing the trajectories on a finite number of basis functions. While raw signals are unlikely to be described by a small number of parameters, this is not the case for smoothed versions of these signals which capture the important patterns. In particular, given a family of smoothed observed trajectories, one may suppose that there exists a basis such that the projection error on a certain number of basis functions of any trajectory is negligible. From now on, the trajectories we consider are assumed to belong to a space spanned by a finite number of basis functions. For the sake of simplicity, we assume in addition that all the components of the trajectories can be decomposed on the same basis but with different dimensions. Extension to different bases is straightforward and does not change our findings but burdens the notation.

Definition 2.5. Let {ϕ k } +∞ k=1 be an orthonormal basis of L 2 [0, T ], R with respect to the inner product

f, g = T 0 f (t) g(t) dt , such that each ϕ k is continuous on [0, T ] and let K := {K d } D
d=1 be a sequence of integers with K := D d=1 K d . We define the space of projected trajectories

Y K (0, T ) ⊂ C [0, T ], R D over [0, T ] as Y K (0, T ) := D d=1 span {ϕ k } K d k=1 .
If there is no risk of confusion, we write Y K := Y K (0, T ) for the sake of readability.

Remark 2.6. From the above definition, any projected trajectory y ∈ Y K is associated with a unique vector c = c

(1) 1 , . . . , c

(1)

K 1 , c (2)
1 , . . . , c

(2)

K 2 , . . . , c (D) 1 , . . . , c (D) K D T ∈ R K defined by c (d) k := y (d) , ϕ k = T 0 y (d) (t) ϕ k (t) dt . (2)
In other words, the vector c is the image of the trajectory y by the projection operator Φ : C [0, T ], R D -→ R K defined by Φy := c, whose restriction Φ| Y K is bijective (as the Cartesian product of bijective operators). In particular, the spaces Y K and R K are isomorphic, i.e. Y K R K .

Regarding the endpoints conditions introduced in definition 2.2, we prove in the following result that satisfying these conditions is equivalent to satisfying a linear system for a projected trajectory.

Proposition 2.7. A trajectory y ∈ Y K belongs to D(y 0 , y T ) if and only if its associated vector c := Φy ∈ R K satisfies the linear system

A(0, T ) c = Γ , (3)
where the matrix A(0, T ) ∈ R 2D×K and the vector Γ ∈ R 2D are defined as follows

A(0, T ) :=           ϕ 1 (0) . . . ϕ K 1 (0) . . . ϕ 1 (0) . . . ϕ K D (0) ϕ 1 (T ) . . . ϕ K 1 (T ) . . . ϕ 1 (T ) . . . ϕ K D (T )          
, Γ := y 0 y T .

Proof.Let y ∈ Y K and let c := Φy ∈ R K . By the definition of the matrix A(0, T ), we have

A(0, T ) c = A(0, T ) c (1) 1 , . . . , c (1) K 1 , c (2) 1 , . . . , c (2) K 2 , . . . , c (D) 1 , . . . , c (D) K D T = K 1 k=1 c (1) k ϕ k (0), . . . , K D k=1 c (D) k ϕ k (0), . . . , K 1 k=1 c (1) k ϕ k (T ), . . . , K D k=1 c (D) k ϕ k (T ) T = y(0) y(T ) .
The conclusion follows directly from the preceding relation.

Reference trajectories modelling

Let us now suppose that we have access to I recorded trajectories y R 1 , . . . , y R I , called reference trajectories, coming from some experiments. We propose here a statistical modelling for these reference trajectories, permitting especially to exhibit some linear properties. This modelling will permit to take advantage of the information contained in these recorded trajectories when deriving optimisation problems in the next subsection.

These trajectories being recorded, they are in particular admissible and we assume that they belong to the space Y K (0, T ). As explained previously they may be interpreted as smoothed versions of recorded signals. In particular each reference trajectory y R i is associated with a unique vector c R i ∈ R K . Moreover we consider each reference trajectory as a noisy observation of a certain admissible and projected trajectory y * . In other words we suppose that there exists a trajectory y * ∈ Y K ∩ A G (y 0 , y T ) associated with a vector

c * ∈ R K satisfying ∀ i = 1, . . . , I c R i = c * + ε i .
The noise ε i is here assumed to be a centered Gaussian whose covariance matrix Σ i is of the form

Σ i = 1 2ω i Σ ,
where Σ ∈ R K×K . It is noteworthy that this matrix will not be known in most of the cases but an estimated covariance matrix can be computed on the basis of the reference vectors. The positive real numbers ω i are here considered as weights so we require I i=1 ω i = 1 ; each ω i plays actually the role of a noise intensity. Further from the hypothesis that the trajectory y * and all the reference trajectories y R i verify the same endpoints conditions, we deduce

A c R i = A c * + A ε i ⇐⇒ A ε i = 0 R 2D ⇐⇒ ε i ∈ ker A ,
for all i = 1, . . . , I (we shorten A(0, T ) in A when the context is clear). Hence the reference vector c * satisfies the following I systems:

       c R i = c * + ε i ε i ∼ N(0 R K , Σ i ) ε i ∈ ker A . ( 4
)
To establish a more explicit system which is equivalent to the preceding one, we require the following preliminary proposition. Here we diagonalise the matrices Σ and A T A by exploiting the fact that the image of the first one is contained in the null space of the other one and vice versa; this is shown in the proof. This property is actually a consequence of the above modelling: the endpoints conditions modelled by A imply linear relations within the components of the vectors, which should be reflected by the covariance matrix Σ. The following result will be helpful to establish proposition 2.10.

Proposition 2.8. We define σ := rank Σ and a := rank A T A. In the setting of system (4), we have σ + a K and there exist an orthogonal matrix V ∈ R K×K and two matrices Λ Σ ∈ R K×K and Λ A ∈ R K×K of the following form:

Λ Σ = Λ Σ,1 0 R σ×(K-σ) 0 R (K-σ)×σ 0 R (K-σ)×(K-σ) , Λ A = 0 R (K-a)×(K-a) 0 R (K-a)×a 0 R a×(K-a) Λ A,2
, where Λ Σ,1 ∈ R σ×σ and Λ A,2 ∈ R a×a are diagonal matrices with positive elements, such that

Σ = V Λ Σ V T , A T A = V Λ A V T .
Proof.The starting point of the proof is to remark that we have

Σ A T A = A T A Σ = 0 R K×K .
(5) Indeed using the hypothesis ε i ∈ ker A for any i = 1, . . . , I gives

Σ A T A = 2ω i Σ i A T A = 2ω i E(ε i ε T i ) A T A = 2ω i E ε i (Aε i ) T A = 0 R K×K ;
similar arguments prove the second equality in (5). First, we can deduce

Im Σ ⊆ ker A T A , (6)
which leads to σ K -a by the rank-nullity theorem. Equalities (5) show also that Σ and A T A are simultaneously diagonalisable (since they commute) so there exists an orthogonal matrix

V ∈ R K×K such that Σ = V Λ Σ V T , A T A = V Λ A V T , (7)
where Λ Σ ∈ R K×K and Λ A ∈ R K×K are diagonal matrices. Permuting if necessary columns of V , we can write the matrix Λ Σ as follows:

Λ Σ = Λ Σ,1 0 R σ×(K-σ) 0 R (K-σ)×σ 0 R (K-σ)×(K-σ) ; ( 8
)
in other words the σ first column vectors of V span the image of Σ. From the inclusion (6), we deduce that these vectors belong to the null space of A T A. Hence the σ first diagonal elements of Λ A are equal to zero and, up to a permutation of the K -σ last column vectors of V , we can write

Λ A = 0 R (K-a)×(K-a) 0 R (K-a)×a 0 R a×(K-a) Λ A,2
, which ends the proof.

Remark 2.9. From equalities (5), we can also deduce

Im A T A ⊆ ker Σ ,
showing that Σ is singular. Consequently the Gaussian noise ε i involved in (4) is degenerate.

A new formulation of system (4) which makes explicit the constrained and unconstrained parts of a vector satisfying this system is given in the following result. This is achieved by using the preceding result which allows to decompose the space R K into three orthogonal subspaces. We prove that the restriction of the noise ε i to the first subspace is a nondegenerate Gaussian, showing that this first subspace corresponds to the unconstrained one. The two other subspaces describe affine relations coming from the endpoints conditions and from implicit relations within the vector components. These implicit relations, which may model for instance natural trends, are expected to be contained in the reference vectors c R i and reflected by the (estimated) covariance matrix Σ. Prior to this, let us write the matrix V ∈ R K×K introduced in proposition 2.8 as follows:

V = V 1 V 2 V 3 , where V 1 ∈ R K×σ , V 2 ∈ R K×K-σ-a and V 3 ∈ R K×a .
We emphasise that the columnvectors of the matrices V 1 and V 3 do not overlap according to the property σ + a K proved in proposition 2.8. In particular the matrix V 2 has to be considered only in the case σ + a < K. Further for any c ∈ R K , we will use the notations

c := V T c , c := V T c ,
for = 1, 2, 3. Finally we consider the singular value decomposition of A coming from the diagonalisation of the symmetric matrix A T A with V :

A = U S A V T ,
where U ∈ R 2D×2D is orthogonal and S A ∈ R 2D×K is a rectangular diagonal matrix of the following form:

S A = 0 R 2D×K-2D S A,2 , (9)
with

S A,2 := Λ A,2 ∈ R 2D×2D .
Proposition 2.10. Suppose that the matrix A is full rank, i.e. a = 2D. Then for any i = 1, . . . , I, system (4) is equivalent to the following one:

               c R i ,1 = c * ,1 + ε i,1 ε i,1 ∼ N 0 R σ , 1 2ω i Λ Σ,1 c * ,2 = V T 2 c R i c * ,3 = S -1 A,2 U T Γ . ( 10
)
Proof.We first prove that system (4) is equivalent to

         c R i = c * + ε i ε i ∼ N 0 R K , 1 2ω i Λ Σ S A c * = U T Γ . ( 11
)
The matrix V being orthogonal, it is non-singular and so we have for all i = 1, . . . , I,

c R i = c * + ε i ⇐⇒ c R i = c * + ε i , and, since Σ i = 1 2ω i Σ = 1 2ω i V Λ Σ V T , we obtain ε i ∼ N(0 R K , Σ i ) ⇐⇒ ε i ∼ N 0 R K , 1 2ω i Λ Σ .
Finally the property ε i ∈ ker A is equivalent to

A c * = Γ ⇐⇒ U S A V T c * = Γ ⇐⇒ S A c * = U T Γ ,
proving that the systems (4) and ( 11) are equivalent. Now the fact that the K -σ last diagonal elements of Λ Σ are zero implies that the components c * ,2 ∈ R K-σ-2D and c * ,3 ∈ R 2D are constant. From the first equality of (11), we have on one side

c R i ,2 = c * ,2 ⇐⇒ V T 2 c R i = c * ,2 ,
for any i = 1, . . . , I. On the other side, combining the last relation of the system (11) with the form of the matrix S A given in (9) permits to obtain

S A c * = U T Γ ⇐⇒ S A,2 c * ,3 = U T Γ ⇐⇒ c * ,3 = S -1 A,2 U T Γ ,
the last equivalence being justified by the hypothesis that the matrix A is full rank (which implies that the diagonal matrix S A,2 is non-singular).

The above decomposition gives us access to non-degenerated density of c R i ,1 given c * ,1 which is later denoted by u( c R i ,1 | c * ,1 ). In next section we will assume a prior distribution on c * ,1 with high density for low values of the cost function F .

A trajectory optimisation problem via a Maximum A Posteriori approach

Before introducing the Bayesian framework, let first recall that we are interested in minimising a certain cost function F : C [0, T ], R D -→ R over the set of projected and admissible trajectories Y K ∩A G (y 0 , y T ). As explained previously, we propose here a methodology leading to a constrained optimisation problem based on the reference trajectories and designed to provide realistic trajectories. Technically speaking, we seek for the mode of a posterior distribution which contains information from the reference trajectories. The aim of this subsection is then to obtain the posterior distribution via Bayes's rule, using in particular the precise modelling of the reference trajectories given in proposition 2.10 and defining an accurate prior distribution with high density for low values of the cost function F .

To do so, we recall firstly that all the trajectories considered here are assumed to belong to the space Y K which is isomorphic to R K . So each trajectory is here described by its associated vector in R K , permitting in particular to define distributions over finitedimensional spaces. We recall also that the reference trajectories are interpreted as noisy observations of a certain y * associated with a c * . According to proposition 2.10, this vector complies with some affine conditions which are described by the following subspace V 1 :

c ∈ V 1 ⇐⇒ V T 2 c = V T 2 c R i V T 3 c = S -1 A,2 U T Γ . ( 12
)
Hence a vector c belonging to V 1 is described only through its component c 1 := V T 1 c. In addition we note that the definition of V 1 does not depend actually on the choice of i since V T 2 c R i has been proved to be constant in proposition 2.10. Further we emphasise that the matrix A is supposed to be full rank in this case and we have V 1 R σ ; we recall that σ is the rank of the covariance matrix Σ.

Let us now define the cost function F over the spaces R K and V 1 . This is necessary to define the prior distribution and to establish our optimisation problem.

Definition 2.11 (Cost functions). Let F : R K -→ R and F : R σ -→ R be the functions defined by

• F (c) := F Φ| -1 Y K c ; • F ( c 1 ) := F Φ| -1 Y K V c T 1 c T R i V 2 Γ T U S -1 A,2 T T .
Remark 2.12. From the preceding definition, we observe that for any y ∈ Y K and its associated vector c ∈ R K , we have

F (c) = F Φ| -1 Y K c = F (y) .
Further for any c ∈ V 1 , we have

F (c) = F Φ| -1 Y K c = F Φ| -1 Y K V c = F Φ| -1 Y K V c T 1 c T R i V 2 Γ T U S -1 A,2 T T = F ( c 1 ) .
We deduce that F is actually the restriction of F to the subspace V 1 .

From now on, the trajectory y * and the associated vector c * will be considered as random variables and will be denoted by y and c. We are interested in the posterior distribution

u( c 1 | c R 1 ,1 , . . . , c R I ,1 ) ,
which depends only on the free component c 1 of c ∈ V 1 , the two other ones c 2 and c 3 being fixed according to (12). We use Bayes's rule to model the posterior via the prior and likelihood distributions, leading to

u( c 1 | c R 1 ,1 , . . . , c R I ,1 ) ∝ u( c R 1 ,1 , . . . , c R I ,1 | c 1 ) u( c 1 ) .
Assuming now that the vectors c R i ,1 are independent gives

u( c R 1 ,1 , . . . , c R I ,1 | c 1 ) u( c 1 ) = I i=1 u( c R i ,1 | c 1 ) u( c 1 ) .
The above likelihood is given by the modelling of the reference trajectories detailed in proposition 2.10. In this case, we have

u( c R i ,1 | c 1 ) ∝ exp -ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 .
The prior distribution is obtained by assuming that the most efficient trajectories (with respect to the cost function) are a priori the most likely ones:

u( c 1 ) ∝ exp -κ -1 F ( c 1 ) , (13)
where κ > 0. Putting everything together and taking the negative of the logarithm gives the following minimisation problem, whose solution is the Maximum A Posteriori estimator:

           c 1 ∈ arg min c 1 ∈R σ F ( c 1 ) + κ I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 c 2 = V T 2 c R i c 3 = S -1 A,2 U T Γ , ( 14
)
where i is arbitrarily chosen in {1, . . . , I}.

Let us now rewrite the above optimisation problem with respect to the variable c = V c ∈ R K in order to make it more interpretable.

Proposition 2.13. The optimisation problem (14) is equivalent to the following one:

c ∈ arg min c∈V 1 F (c) + κ I i=1 ω i c -c R i T Σ † c -c R i , (15)
where Σ † ∈ R K×K denotes the pseudoinverse of the matrix Σ.

Proof.From (8), we deduce

I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 = I i=1 ω i c -c R i T Λ † Σ c -c R i = I i=1 ω i c -c R i T V Λ † Σ V T c -c R i = I i=1 ω i c -c R i T Σ † c -c R i .
And from the proof of proposition 2.10, we have

A c = Γ ⇐⇒ c 3 = S -1 A,2 U T Γ , proving that c ∈ V 1 .
To conclude, let us comment on this optimisation problem.

1. To interpret the optimisation problem (15) (or equivalently ( 14)) from a geometric point of view, let us consider the following new problem:

min c 1 ∈R σ F ( c 1 ) s.t. I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 κ (16)
where λ 0. Here we suppose that F is strictly convex and that the problem (16) has a solution (which is then unique). By Slater's theorem (Boyd and Vandenberghe, 2004, Subsec. 5.2.3), the strong duality holds for the problem (16). It can then be proved that there exists a certain λ 0 such that the solution of ( 16) is the minimiser of the strictly convex function

c 1 -→ F ( c 1 ) + λ I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 ,
which is actually the objective function of the optimisation problem ( 14) for κ = λ . Hence the problem (14) minimises the cost F in a ball centered on the weighted average of the reference trajectories. In particular if the reference trajectories are close to an optimal one with respect to F then one could expect the solution of ( 14) to be equal to this optimal trajectory.

2. Further the optimisation problem (15) takes into account the endpoints conditions through the subspace V 1 but not the additional constraints. However as explained in the preceding point, the solution is close to realistic trajectories and so is likely to comply with the additional constraints for a well-chosen parameter κ > 0. We refer to Sec. 2.6 for more details on an iterative method for the tuning of κ. In particular, a right choice for this parameter is expected to provide an optimised trajectory with a realistic behaviour. This is for instance illustrated in Sec. 4.

3. Taking into account the linear information from the available data through the covariance matrix Σ allows to restrict the search to the subspace V 1 describing these relations. This is of particular interest when implicit relations (modelled by the sub-matrix V 2 ) are revealed by the estimation of Σ on the basis of the reference trajectories; in this case, these implicit relations may not be known by the expert.

4. The optimisation problem (15) has linear constraints and a quadratic penalised term.

For instance, if the cost function F is a convex function then we obtain a convex problem for which efficient algorithms exist.

Quadratic cost for a convex optimisation problem

In this short subsection, we focus on a particular case where the cost function F is defined as the integral of an instantaneous quadratic cost function, i.e.

∀ y ∈ C [0, T ], R D F (y) = T 0 f (y(t)) dt , (17)
where f : R D -→ R is quadratic. Even though such a setting may appear to be restrictive, we emphasise that quadratic models may lead to highly accurate approximations of variables, as it is illustrated in Sec. 4. For a quadratic instantaneous cost, the associated function F : R K -→ R can be proved to be quadratic as well and can be explicitly computed. In the following result, we provide a quadratic optimisation problem equivalent to (15).

Proposition 2.14. Suppose that the cost function F is of the form (17) with f quadratic.

Then the optimisation problem (15) is equivalent to the following one:

c ∈ arg min c∈V 1 c T Q + κ Σ † c + w -2κ I i=1 ω i Σ † c R i T c , (18)
where Q ∈ R K×K and w ∈ R K can be explicitly computed from f .

Proof.We defer the proof to the supplementary material.

In particular this permits to derive sufficient conditions on the parameter κ > 0 so that the optimisation problem is proved to be equivalent to a quadratic program (Boyd and Vandenberghe, 2004, Sec. 4.4), namely the objective function is convex quadratic together with affine constraints. In practice, this allows to make use of efficient optimisation libraries to solve numerically (18).

Iterative process to comply with additional constraints

As explained in Sec. 2.4, the trajectory optimisation problem ( 15) is constrained by the endpoints conditions and by implicit linear relations revealed by the reference trajectories. Nevertheless the additional constraints introduced in definition 2.3 are not taken into account in this problem. In practice such constraints assure that natural or user-defined features are verified and so a trajectory which does not comply with these constraints may be considered as unrealistic.

Our aim is then to assure that the trajectory y = Φ| -1 Y K c , where c ∈ V 1 is the solution of the optimisation problem (15), verifies the additional constraints, i.e. belongs to the set G. A first solution would be to add the constraint Φ| -1 Y K c ∈ G in the optimisation problem (15). However depending on the nature of the constraints functions g , this may lead to nonlinear constraints which could be costly from a numerical point of view. The solution we propose consists rather in exploiting the degree of freedom coming from the parameter κ > 0 appearing in the problem (15).

First of all, let us factorise the problem (15) by κ to obtain the following new one for the sake of presentation:

c ∈ arg min c∈V 1 ν F (c) + I i=1 ω i c -c R i T Σ † c -c R i , (19)
where ν := κ -1 . On one hand, we observe that the solution of the optimisation problem (19) for the limit case ν = 0 is given by I i=1 ω i c R i which is the average of the reference vectors. In this case, one may expect that the associated average trajectory complies with the constraints but is unlikely to optimise the cost function F . On the other hand, for very large ν > 0, the second term of the objective function in ( 19) can be considered as negligible compared to the first one. In this case, the cost of the solution will surely be smaller than the costs of the reference trajectories but no guarantee regarding the additional constraints can be established in a general setting.

Given these observations, the task is then to find an appropriate value ν > 0 in order to reach a trade-off between optimising and remaining close to the reference trajectories to comply with the additional constraints. Many methods can be developed to find such a ν and, among those based on iterative processes, linear or binary search algorithms can be considered. In this case, one has to set firstly a maximal value ν max so that the solution of ( 19) with ν max is unlikely to satisfy the constraints and to perform secondly the search over the interval (0, ν max ). Since the solution for ν = 0 is assumed to be admissible, we expect that the binary search will find a ν > 0 leading to an optimised trajectory belonging to G.

Confidence bounds on the integrated cost

In practice the cost function F considered is an estimation of the true cost F , a random variable which cannot be fully predicted based on y. If the distribution F (y) would be known it would be possible to deduce confidence bound on F . This is for instance possible by considering multivariate functional regression [START_REF] Ramsay | Parameter estimation for differential equations: a generalized smoothing approach[END_REF].

The simplest case from the estimation point of view is to consider that F is the integral of some instantaneous consumption function f as in Sec. 2.5, and to estimate the parameters of the standard multivariate regression

f (y(t)) = f (y(t)) + ε(t),
where the random noise ε(t) is assumed to follow a centered Gaussian distribution with variance σ. In this case F can be expressed as the integral of a stochastic process

F (y) := T 0 f (y(t)) dt = F (y) + T 0 ε(t) dt .
then assuming that (ε(t)) t∈[0,T ] independent we obtain

T 0 ε(t) dt ∼ N(0, T σ 2 ).
Thus F (y) follows a Gaussian distribution centered on F (y) and with variance equals to T σ 2 . This makes it possible to compute confidence bounds on F (y). For a confidence level 1 -u, u ∈ [0, 1], a confidence interval for F (y) is obtained as

CI 1-u (F (y)) = F (y) ± ζ 1-u 2 √ T σ,
where ζ 1-u 2 is the quantile of order 1 -u 2 of the standard Gaussian distribution. The assumption that f and σ 2 are known is relevant since they are estimated based on a huge amount of training data. The assumption of white Gaussian noise can be seen as unrealistic, however it appears to be the only route to explicit calculus. A more complex strategy could be derived using Gaussian processes, which is beyond the scope of this paper.

The Python library Pyrotor

The above optimisation methodology is aimed at being used in a wide range of applications, from path planning for industrial robots [START_REF] Chettibi | Minimum cost trajectory planning for industrial robots[END_REF] to fuel-efficient aircraft trajectories [START_REF] Dewez | From industry-wide parameters to aircraft-centric on-flight inference: improving aeronautics performance prediction with machine learning[END_REF][START_REF] Rommel | Gaussian mixture penalty for trajectory optimization problems[END_REF]. We therefore contribute a generic Python library PyRotor (standing for Python Route trajectory optimiser) which is intended to a large audience. In particular, this library has been used to obtain the numerical results given in the two following sections.

When using the PyRotor library, the practitioner has to define the endpoints conditions as a dictionary, the additional constraints in a list of functions, the name of the basis and the dimension for each variable. The current version of the library covers only the case of Sec. 2.5, that is to say the cost is given by a quadratic instantaneous function. This permits to make use of proposition 2.14 in which a quadratic objective function is given. We mention that future releases of PyRotor are intended to cover more general cost functions. The value of the parameter ν max in (19) can also be manually set depending on the application. The Legendre basis is currently the only basis implemented in the first version of PyRotor (via the legendre module from NumPy package [START_REF] Harris | Array programming with NumPy[END_REF]) but future developments including other general bases are planned. Further the user indicates a path to a directory containing the data, each reference trajectory being contained in a csv file. The covariance matrix Σ is here estimated by using the sklearn.covariance package from the Python library scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Two optimisation solvers are proposed: the generic solver minimize(method='trust-constr') [START_REF] Conn | Trust Region Methods[END_REF] from SciPy software [START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] and the quadratic programming solver from CVXOPT software [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF]. The latter is intended to speed up the execution in case of convex quadratic objective function. Once the arguments are given by the user, a class is created and the optimisation is performed by executing a method from this class. At the end, the optimised trajectory is provided in a dataframe: at each time, the position of the trajectory is given together with the value of f . The total cost is also computed and a quantitative comparison in terms of savings with the reference trajectories can be also displayed.

The open source PyRotor library is developed on GitHub and welcomes contributions from its users: we favour a community-based development to foster the diffusion of our work towards practitioners. PyRotor is intended to be PEP8 compliant and purposely rely on high standard coding practices. The continuous development platform Travis is used to certify the latest builds of the library. Finally we provide Jupyter notebooks for examples on how to use PyRotor along with online documentation. PyRotor is available at https://github.com/bguedj/pyrotor.

Application 1: trajectory optimisation for fuel-efficient aircrafts

In this section, we consider the aeronautic problem of reducing the total fuel consumption of an aircraft during the climb phase. This example illustrates the key role played by the reference trajectories since we are able to obtain realistic and optimised trajectories thanks to a simple modelling involving few constraints.

Modelling

Here the trajectories are supposed to be in a vertical plane and are defined by the altitude h, the Mach number M and the engines rotational speed N1 (expressed as a percentage of a maximal value). Hence a trajectory y in this setting is a continuous R 3 -valued map defined on [0, T ], where T is a maximal climb duration fixed by the user. Hence we have

∀ t ∈ [0, T ] y(t) := h(t), M(t), N1(t) .
The quantity to minimise is the total fuel consumption TFC : Regarding the endpoints conditions, we require the trajectory to start at the altitude h 0 with Mach number M 0 and to end at the altitude h T with Mach number M T . In particular, the reference trajectories we use have to verify these conditions. We consider also additional constraints which are conventional in the aeronautic setting:

C [0, T ], R 3 -→ R +
• The rate of climb, i.e. the time-derivative of the altitude, has to be upper bounded by a given maximal value γ max during the whole climb;

• The Mach number should not exceed a certain value called the maximum operational Mach (MMO).

The final time of the climb is given by T ∈ [0, T ] which is the first time where the aircraft reaches h T with Mach number M T . Finally we mention that the fuel flow model FF is here estimated. To do so, we exploit the reference trajectories which contain recorded altitude, Mach number, engines power and fuel flow for each second of the flight. Having access to these data, we are in position to fit a statistical model. Following the numerical results in [START_REF] Dewez | From industry-wide parameters to aircraft-centric on-flight inference: improving aeronautics performance prediction with machine learning[END_REF] which show that polynomials can accurately model aeronautic variables, we consider a polynomial model of degree 2 for the fuel flow. In particular the requirements for the cost function in the current version of PyRotor are fulfilled. The prediction accuracy of the resulting estimated model is assessed in the following subsection.

Numerical results

We present now numerical results based on real flight data for the above aeronautic problem. Here we have access to 2,162 recorded short and medium-haul flights performed by the same narrow-body airliner type, provided by a partner airline. In particular they can not be publicly released for commercial reasons. The data is here recorded by the Quick Access Recorder (QAR).

Before considering the optimisation setting, we estimate a fuel flow model specific to the climb phase and to the considered airliner type. To do so we extract the signals of the four variables of interest (altitude, Mach number, engines rotational speed and fuel flow) and keep the observations from the take-off to the beginning of the cruise without level-off phases. Smoothing splines are then applied to the raw signals to remove the noise. We sample each 5 seconds to reduce the data set size without impacting strongly the accuracy of the resulting models. At the end, we obtain 494,039 observations which are randomly split into training and test sets to fit a polynomial model of degree 2 using the scikit-learn library. The RMSE and MAPE values of this model on the test set are respectively equal to 3.64 × 10 -2 kg.s -1 and 1.73%.

Regarding the optimisation, we are interested in climb phases from 3,000 ft to 38,000 ft. We mention that we remove lower altitudes because operational procedures constraint heavily the trajectory during the very beginning of the climb. Further the initial and final Mach numbers are required to be equal to 0.3 and 0.78. It is noteworthy that the optimisation solvers used in PyRotor allow linear inequality conditions, permitting to slightly relax the endpoints conditions. Here we tolerate an error of 100 ft for the altitude and an error of 0.01 for the Mach number. The initial and final N1 values are let unconstrained. Finally the MMO and γ max are respectively set to 0.82 and 3,600 ft.min -1 .

The reference trajectories are given by 48 recorded flights which satisfy the above climb endpoints conditions among the 2,162 available ones. All these selected flights are used to estimate the covariance matrix involved in the optimisation problem. On the other hand, we use only the 5 most fuel-efficient flights in the objective function to focus on a domain containing the most efficient recorded flights. Further the maximal duration T is here fixed to the duration of the longest climb among the 5 most fuel-efficient ones we use.

Legendre polynomials are used as the functional basis spanning the space in which lies the trajectories. Since we consider narrow-body airliners, polynomials are expected to be relevant to describe the slow variations of such aircrafts. Here the dimensions associated with the altitude, the Mach number and the engines power are given respectively by 4, 10 and 6. The reference vectors c R i are then computed using the formula (2). At the end, we amount to solving a constrained optimisation problem in a space of dimension 20.

We are then in position to apply the optimisation method developed in Sec. 2 using the PyRotor library. First of all a relevant value for ν max > 0 has to be fixed. In order to propose a realistic optimised climb, we choose a ν max relatively small so that the optimised climb remains close to the reference ones. In particular, the quadratic objective function in (19) turns out to be convex for all ν ∈ [0, ν max ] permitting to use the quadratic programming solver from CVXOPT software imported in PyRotor. The preprocessing of the reference trajectories and the optimisation steps have been executed 100 times using PyRotor on an Intel Core i7 6 cores running at 2.2 GHz. The mean of the execution time for both steps is equal to 3.76 s with standard deviation 0.11 s, illustrating that the library is time-efficient in this setting.

A plot of the optimised trajectory obtained using PyRotor is given in Fig. 1. We observe that the optimised trajectory seeks to reach the maximum altitude in the minimum amount of time; this is in accordance with the existing literature (see for instance Codina Table 1: Statistical description of the fuel savings of the optimised trajectory -The savings are compared with the 48 recorded flights satisfying the present endpoints and the total consumption of the optimised trajectory is estimated by using the statistical model for the fuel flow -Q 1 , Q 2 and Q 3 refer to the first, second and third quartiles. and Menéndez (2014) and references therein). In particular, the duration T is equal to 1,033 seconds which is actually slightly shorter than the reference durations. We note also that the optimised Mach number shares a very similar pattern with the references.

Mean Std Min Q 1 Q 2 Q 3 Max Fuel
On the other hand, the optimised engines rotational speed tends to slowly decrease until the cruise regime before reaching the top of climb. This is not the case for the reference engines speed which falls to the cruise regime just after reaching the final altitude. Most of the savings seem to be achieved in these last moments of the climb. At last but not least, the optimised trajectory presents a realistic pattern inherited from the reference trajectories.

For a quantitative comparison, we refer to Table 1 which provides statistical information on the fuel savings. The mean savings 16.54% together with the fact that the optimised trajectory verifies the additional constraints show that these first results are promising, motivating further studies. For instance one could model environmental conditions or take into account Air Traffic Control constraints for more realistic modellings.

Application 2: trajectory optimisation to maximise work of a force field

Here we consider the following generic example: given a moving point in a force field, find a trajectory starting and ending at two different given points which maximises the work of the force along the trajectory while minimising the travelled distance. For instance, this corresponds to a very simple modelling of a sailing boat which seeks to increase the power of the wind at each time, i.e. maximising the wind work, without travelling a too large distance. This second example demonstrates that our generic optimisation approach is flexible enough to take into account derivatives of trajectories and hence to cover dynamics settings.

Modelling

To model this problem, we suppose without loss of generality that the trajectories are defined on the (time-)interval [0, 1] and we let V : R D -→ R D denote a vector field. Furthermore the trajectories are assumed here to be continuously differentiable, i.e. they belong to

C 1 [0, 1], R D . The work of V along a trajectory y ∈ C 1 [0, 1], R D is defined as W (y, ẏ) := 1 0 V y(t) T ẏ(t) dt ;
here ẏ denotes the derivative of y with respect to the independent variable t. Moreover using Hamilton's principle in Lagrangian mechanics, it can be shown that the trajectory with constant velocity (i.e. a straight line travelled at constant speed) is the minimum of the following functional,

J( ẏ) = 1 0 ẏ(t) 2 2 dt ,
where the starting and ending points of y are fixed and different. This functional can be then used to control the travelled distance. It follows that minimising the cost function

F α (y, ẏ) := αJ( ẏ) -W (y, ẏ) = 1 0 α ẏ(t) 2 2 -V y(t) T ẏ(t) dt ,
where α 0 is arbitrarily chosen, is expected to lead to an optimised trajectory reflecting a trade-off between maximising the work and minimising the distance. Further we require the trajectory to stay in the hypercube [0, 1] D and to start and to end respectively at y 0 ∈ [0, 1] D and y 1 ∈ [0, 1] D . Now we remark that the above cost function involves the (time-)derivative ẏ. So one has to derive a formula permitting to compute the derivative of any trajectory y = Φ| -1 Y K c ∈ Y K from its associated vector c ∈ R K , especially to compute F (c). For instance, this can be easily achieved by assuming that each element of the functional basis is continuously differentiable. Indeed we can differentiate in this case any y ∈ Y K :

∀ d = 1, . . . , D ẏ(d) = K d k=1 c (d) k φk = d dt Φ| -1 Y K c (d)
.

We deduce then the following formula for F (c) in the present setting:

F (c) := F α Φ| -1 Y K c, d dt Φ| -1 Y K c .
Here the vector c contains information on both position and velocity, permitting especially to keep the problem dimension unchanged. To finish, let us remark that it is possible to make the above formula for F explicit with respect to c in certain settings. For instance it is possible to derive an explicit quadratic formula for F (c) when the integrand defining F α is quadratic with respect to y(t) and ẏ(t); this formula is implemented in PyRotor and the arguments to obtain it are similar to those proving proposition 2.14.

Numerical results

Numerical results based on randomly generated data for the above physical application are presented in this section.

First of all we consider trajectories with two components y (1) and y (2) lying in the square [0, 1] 2 for the sake of simplicity. We set the starting and ending points as follows: y (1) (0) = 0.111 , y (2) (0) = 0.926 , y (1) (1) = 0.912 , y (2) (1) = 0.211 with a tolerated error 1 × 10 -4 , and the vector field V : R 2 -→ R 2 is here defined by V x (1) , x (2) = 0, x (1) T .

Given the above endpoints and the vector field, we observe that the force modelled by V will be in average a resistance force to the motion. Indeed the force is oriented toward the top of the square while the moving point has to go downward. Further let us note that the Here the reference trajectories are obtained through a random generation process. To do so, we define an arbitrarily trajectory y R verifying the endpoints conditions and we compute its associated vector c R ; Legendre polynomials are once again used and the dimensions of y (1) and y (2) are here set to 4 and 6. Let us note that y R is designed in such a way that it has a relevant pattern but not the optimal one. Then we construct a set of reference trajectories by adding centered Gaussian noises to c R . It is noteworthy that the noise is generated in such a way that it belongs to the null space of the matrix A describing the endpoints conditions; the resulting noised trajectories satisfy then these conditions. Further the trajectories which go out of the square [0, 1] 2 are not kept. At the end, we get 122 generated reference trajectories assumed to be realistic in this setting, each of them containing 81 time observations. Among these reference trajectories, we use the 10 most efficient ones with respect to the cost F α .

In the present example, we set a ν max relatively large to explore a large domain around the reference trajectories. In this case, the objective function of the optimisation problem (19) may be not convex even if it is still quadratic. So we make use of the generic optimisation solver minimize(method='trust-constr') imported in PyRotor. Regarding the execution time, we have randomly and uniformly generated 100 values in the interval [0, 10] for the parameter α and executed PyRotor for each of them. The mean of PyRotor execution time is 0.44 s with standard deviation 0.03 s on an Intel Core i7 6 cores running at 2.2 GHz.

In Fig. 2, we plot 4 optimised trajectories associated with different values of α: 0, 0.35, 1 and 10. As expected the trajectory associated with the largest value of α gives the most straight trajectory while the most curvy one is associated with α = 0. In particular, the latter tends to move to the left at the beginning where the force V is the smallest before going to the ending point in a quasi-straightforward way so that the force is perpendicular to the motion. This example illustrates especially that our optimisation approach may lead to optimised trajectories which differ from the reference ones to reduce more the cost.

A quantitative comparison in terms of work gains for different values of α is provided in Table 2. The results confirm the above observations on the curves and show that a right value for α has to be fixed depending on the setting.

Conclusion / Discussion

We have proposed an approach for data-driven trajectories optimisation without involving dynamical system. The approach can work based on a known cost function or a cost function learnt from the data. The modelling of the trajectories allows to take into account explicit and implicit linear constraints on the coefficients in the optimisation problem. Contrary to full optimisation approaches, our method finds a trade-off between high density constraints-compliant solutions and fully optimised solutions through the tuning of the regularisation. In the aeronautic framework our approach leads to promising fuel-efficient trajectories. Our approach is generic enough to be applied to other physical settings such as the motion of a moving point in a force field (such as a sailing boat). Some perspectives of this work are first to further exploit the flexibility of Bayesian setting, by not only searching for the mode of the posterior distribution but also sampling by means of MCMC algorithms. A second perspective would be to consider a clustering of reference trajectories, and apply our strategy on each cluster, then particularise the optimal trajectory depending on the cluster. Last but not least, we aim at adapting our approach where some component of the trajectory would be categorical variables: this would be particularly useful for decision making processes in various disciplines.

cial Intelligence (Grant: 815914). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Data that support the findings of this study have been gathered from the Quick Access Recorder (QAR) of aircrafts operated by a private airline, hence can not be released publicly for commercial reasons. The authors are grateful to Baptiste Gregorutti and Pierre Jouniaux for fruitful discussions about the modelling of the problem and the validation of the results for the aeronautic application.

A. Appendix

Here we focus on the case where the cost function

F : C [0, T ], R D -→ R D is of the following form F (y) = T 0 y(t) T Q y(t) + w T y(t) + r dt , (20)
where Q ∈ R D×D is symmetric, w ∈ R D and r ∈ R. In this setting, we provide explicit formulas for the costs F : R K -→ R and F : R σ -→ R defined in Section 2.4. A sufficient condition on the parameter κ > 0 so that the optimisation problem

c ∈ arg min c∈V 1 F (c) + κ I i=1 ω i c -c R i T Σ † c -c R i , (21)
is a quadratic program in the present setting is then derived. From Section 2.4, we recall that the preceding optimisation problem is equivalent to

           c 1 ∈ arg min c 1 ∈R σ F ( c 1 ) + κ I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 c 2 = V T 2 c R i c 3 = S -1 A,2 U T Γ . ( 22
)
Lemma A.1. Suppose that the cost function F is of the form (20). Then the costs F and F are quadratic and explicit formulas are given in (24) and (26).

Proof.Let c ∈ R K and let y := Φ| -1 Y K c ∈ Y K be its associated trajectory, which can be represented as follows:

∀ d ∈ {1, . . . , D} y (d) = K d k=1 c (d) k ϕ k .
We also remark that each component of the vector

c = c (1) 1 , . . . , c (1)
K 1 , c (2) 1 , . . . , c (2)
K 2 , . . . , c (D) 1 , . . . , c (
D) K D
T can be simply described by a single parameter so that we can write c = (c 1 , c 2 , . . . , c K ) T .

• Computation of F :

We first insert the preceding representation of y into the above quadratic integrand to obtain:

y(t) T Q y(t) + w T y(t) + r = D d 1 ,d 2 =1 K d 1 k 1 =0 K d 2 k 2 =0 Q d 1 d 2 c (d 1 ) k 1 c (d 2 ) k 2 ϕ k 1 (t)ϕ k 2 (t) + D d=1 K d k=0 w d c (d) k ϕ k (t) + r , (23)
for all t ∈ [0, T ]. The next step of the proof consists in changing the indices of the above sums. To do so, let us define the matrix Q ∈ R K×K and the vector w ∈ R K as

Q :=    Q 11 J K 1 ,K 1 . . . Q 1D J K 1 ,K D . . . . . . Q D1 J K D ,K 1 . . . Q DD J K D ,K D    , w := w 1 J 1,K 1 . . . w D J 1,K D T ,
where J m,n is the all-ones matrix of size m × n. We also introduce the map ϕ ∈

C [0, T ], R K as ϕ(t) := ϕ 1 (t), . . . , ϕ K 1 (t), ϕ 1 (t), . . . , ϕ K 2 (t), . . . , ϕ 1 (t), . . . , ϕ K D (t) T , for all t ∈ [0, T ],
where the ϕ k are the functional basis elements. We are now in position to change the indices in the sums appearing in (23):

D d 1 ,d 2 =1 K d 1 k 1 =0 K d 2 k 2 =0 Q d 1 d 2 c (d 1 ) k 1 c (d 2 ) k 2 ϕ k 1 (t)ϕ k 2 (t) + D d=1 K d k=0 w d c (d) k ϕ k (t) + r = K k 1 ,k 2 =1 Q k 1 k 2 c k 1 c k 2 ϕ k 1 (t)ϕ k 2 (t) + K k=1 w k c k ϕ k (t) + r ,
where we have used the above rewriting of the vector c. Integrating finally over [0, T ] gives F (c) = T 0 y(t) T Q y(t) + w T y(t) + r dt

= K k 1 ,k 2 =1 Q k 1 k 2 T 0 ϕ k 1 (t) ϕ k 2 (t) dt c k 1 c k 2 + K k=1 w k T 0 ϕ k (t) dt c k + rT = K k 1 ,k 2 =1 Qk 1 k 2 c k 1 c k 2 + K k=1 wk c k + rT = c T Qc + wT c + rT , (24)
where we have defined

Qk 1 k 2 := Q k 1 k 2 T 0 ϕ k 1 (t) ϕ k 2 (t) dt , wk := w k T 0 ϕ k (t) dt . (25)
• Computation of F :

By the definition of F given in Section 2.4, we have

F ( c 1 ) = F V c T 1 c T 2,3 T ,
where V has been introduced in Section 2.3 and c 2,3 ∈ R K-σ is defined as follows:

c 2,3 := V T 2 c R i S -1
A,2 U T Γ , here the index i is arbitrarily chosen in {1, . . . , I} since the vector V T 2 c R i has been proved to be independent from i. We introduce now the matrix QV := V T QV and the vector wV := V T w which can be decomposed as follows: , where QV,11 ∈ R σ×σ , QV,12 ∈ R σ×(K-σ) , QV,21 ∈ R (K-σ)×σ , QV,22 ∈ R (K-σ)×(K-σ) , wV,1 ∈ R σ and wV,2 ∈ R K-σ . Given this and the preceding point, we obtain Rearranging the preceding terms and using the fact that QV is symmetric gives

F ( c 1 ) = c T 1 c T 2,3 V T QV c T 1 c T 2,3 T + V T w T c T
F ( c 1 ) = c T 1 Q c 1 + w T c 1 + r , (26)
where

• Q := QV,11 ; (27)
• w := 2 QV,12 c 2,3 + w V,1 ;

• r := c T 2,3 QV,22 c 2,3 + w T V,2 c 2,3 + rT .

(29)

The optimisation problem ( 22) is then equivalent to the following one in the present quadratic setting:

           c 1 ∈ arg min c 1 ∈R σ c T 1 Q c 1 + w T c 1 + κ I i=1 ω i c 1 -c R i ,1 T Λ -1 Σ,1 c 1 -c R i ,1 c 2 = V T 2 c R i c 3 = S -1 A,2 U T Γ . ( 30
)
In the following result, we provide a sufficient condition on the parameter κ > 0 so that the problem (30) is a quadratic program. The proof uses the fact that the symmetric matrix associated with the quadratic objective function is now explicit and given by the sum of two matrices. A perturbation result for matrices is then applied to obtain a bound for κ assuring that the symmetric matrix is positive semidefinite. Proof.We first note that all the eigenvalues of the matrix Σ are non-negative (because Σ is a covariance matrix) and that λ σ+1 = . . . = λ K = 0 (because rank Σ = σ). In particular, the eigenvalue λ 1 is positive.

Standard calculations show that the symmetric matrix associated with the quadratic objective function of the problem (30) is given by M (κ) := Q + κ Λ -1 Σ,1 ∈ R σ×σ . Let µ 1 (κ) µ 2 (κ) • • • µ σ (κ) denote the eigenvalues of M (κ). The goal is then to find a sufficient condition on κ > 0 so that µ σ (κ) is non-negative to assure that M is positive semidefinite. Since M (κ) can be interpreted as a perturbed version of Q, we can apply Weyl's inequality (see for instance [START_REF] Wang | The converse of weyl's eigenvalue inequality[END_REF]) which states µ σ (κ) ρ σ + κ λ 1 .

Then choosing κ such that κ -ρ σ λ 1 implies that µ σ (κ) 0, leading to the result.

For the sake of completeness, we finish by rewriting the problem (30) as a quadratic optimisation problem in V 1 ⊂ R K .

Proposition A.3. Suppose that the cost function F is of the form (20). Then the optimisation problem (21) is equivalent to the following one:

c ∈ arg min c∈V 1 c T Q + κ Σ † c + w -2κ I i=1 ω i Σ † c R i T c .
Proof.It is sufficient to show that the two following objective functions g 1 , g 2 : R K -→ R have the same minima:

• g 1 (c) := F (c) + κ I i=1 ω i c -c R i T Σ † c -c R i ; • g 2 (c) := c T Q + κ Σ † c + w -2κ I i=1 ω i Σ † c R i T c .
Firstly we have by standard calculations,

I i=1 ω i c -c R i T Σ † c -c R i = I i=1 ω i c T Σ † c -2 I i=1 ω i c T R i Σ † c + I i=1 ω i c T R i Σ † c R i = c T Σ † c -2 I i=1 ω i Σ † c R i T c + I i=1 ω i c T R i Σ † c R i ,
for any c ∈ R K , where we have used I i=1 ω i = 1. Combining now this equality with Lemma A.1 implies

g 1 (c) = c T Qc + wT c + rT + κ   c T Σ † c -2 I i=1 ω i Σ † c R i T c + I i=1 ω i c T R i Σ † c R i   = c T Q + κ Σ † c + w -2κ I i=1 ω i Σ † c R i T c + κ I i=1 ω i c T R i Σ † c R i + rT = g 2 (c) + κ I i=1 ω i c T R i Σ † c R i + rT .
Since the two last terms of the last right-hand side do not depend on c, we deduce that the objective functions g 1 and g 2 have the same minima.

  which is defined via the fuel flow FF : R 3 -→ R + as follows 1 :

Figure 1 :

 1 Figure 1: Optimised and reference altitudes, Mach numbers and engines rotational speeds -The optimised trajectory is represented by the blue curves.

Figure 2 :

 2 Figure 2: Optimised trajectories in the square [0, 1] 2 for α ∈ {0, 0.35, 1, 10} -Optimised and reference trajectories are respectively given by plain and dotted curves -Coloured dots indicate the power value of the force at different points of the optimised trajectories and the bar shows the scale -Red arrows represent the pattern of the vector field V .

  11 c 1 + c T 1 QV,12 c 2,3 + c T 2,3 QV,21 c 1 + c T 2,3 QV,22 c 2,3 + w T V,1 c 1 + w T V,2 c 2,3 + rT .

  λ K be respectively the eigenvalues of the symmetric matrices Q and Σ. If κ -ρ σ λ 1 then the optimisation problem (30) is a quadratic program.

Table 2 :

 2 Statistical description of the work gains in percentage for α ∈ {0, 0.35, 1, 10} -The values have been computed by using the 122 available reference trajectories -Negative percentages indicate that no work gains have been obtained -Q 1 , Q 2 and Q 3 refer to the first, second and third quartiles.

		Mean	Std	Min	Q 1	Q 2	Q 3	Max
	α = 0	73.43	2.36	68.63	71.90	73.25	74.67	80.69
	α = 0.35	45.88	4.81	36.09	42.75	45.49	48.39	60.66
	α = 1	-6.12	9.43	-25.31	-12.26	-6.88	-1.20	22.87
	α = 10	-34.54	11.96	-58.87	-42.32	-35.50	-28.30	2.22
	integrand of the cost function F α in the present setting is actually quadratic with respect
	to y(t) and ẏ(t), so that an explicit quadratic formula for F (c) implemented in PyRotor
	is available.

In the notation of Sec.

2.5, FF and TFC play respectively the role of f and F .

Acknowledgements

This work was supported by the European research program Clean Sky 2 under the project PERF-AI: Enhance Aircraft Performance and Optimisation through utilisation of Artifi-


Diego Pereira Botelho

Douglas Martins Araujo

Gérard Meunier

Yves Marechal

Jean-Louis Coulomb

Brahim Ramdane

Olivier Chadebec

Hybrid Natural Element Method -Boundary Element Method for Unbounded Problems

Keywords: boundary element method, coupled numerical methods, natural element method

This paper presents a hybrid approach coupling natural element method (NEM) and boundary element method (BEM) for the treatment of unbounded problems. The goal is to combine the accuracy of the NEM and the ability of the BEM in modeling linear and deformable domains without a mesh. The performance of the proposed scheme is compared to a FEM-BEM coupling in terms of accuracy. Results show that for a given number of degrees of freedom the proposed approach (NEM-BEM) is able to provide more accurate solutions.

I. INTRODUCTION

he Natural Element Method (NEM) is a numerical method based on the discretization by Voronoi diagram and on the concept of natural neighbors [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. This method presents some characteristics of meshfree methods [START_REF] Belytschko | Meshfree Methods[END_REF], like the highly accurate approximations and the shape functions with quasispherical influence zones. However, unlike most of the meshless methods, the NEM shape functions are interpolants, allowing the easy handling of the interface between materials and the direct imposition of essential boundary conditions. In terms of performance, the NEM has been proved to be more efficient (lower computational cost for a given accuracy) than the 1 st order finite element method (FEM) [START_REF] Maréchal | Computational Performances of Natural Element and Finite Element Methods[END_REF].

In this work, in order to model unbounded regions, we use an integral method known as the Boundary Element Method (BEM) [START_REF] Brebbia | The Boundary Element Method For Engineers[END_REF]. Using this approach, the degrees of freedom are located on the interface between the domains (bounded and unbounded).

By combining the accuracy of the natural element method and the ability of the boundary element method in modeling linear and deformable domains without a mesh, the quality of interpolation on the magnetic material and of the field reconstitution in the unbounded region will be improved.

The developed approach is validated in the case of two 2D magnetostatic problems. The results are compared with analytical solution in the first case and with the classical finite element method for the second one.

II. MAGNETOSTATIC PROBLEM

The 2D magnetostatic problem is defined by the Maxwell's equations:

j h = × ∇ , (1) 0 = ⋅ ∇ b . ( 2
)
The representation of the modeled domains is showed in Fig. 1. Ω m is the magnetic region (which can contain currents and magnets), Ω 0 is the unbounded region (may optionally contain non-meshed coils) and Γ m is the interface between magnetic and unbounded domains. In order to introduce the non-meshed coil a reduced magnetic scalar potential φ m is defined, as follow:

0 h h h + = red , ( 3
)
m red ϕ -∇ = h , (4)
where h 0 is the source field generated by the non-meshed coil and h red is the reduced magnetic field.

III. NATURAL ELEMENT METHOD

The natural element method uses the concept of natural neighbors. It is based on the construction of Voronoï diagram on a cloud of nodes. For a set of nodes N = {n 1 , n 2 , n 3 … n N } distributed in a given 2-D domain, the Voronoï diagram is a subdivision of the domain into cells, where each cell C i associated to node n i is such that any point in C i is closer to node n i than to any other node n j for i ≠ j. In mathematical terms, a Voronoï cell is defined as:

( ) ( ) { } i j , , d ,d C n i ≠ ∀ < ℜ ∈ = j i x x x x x , , (5)
where d is the distance between two points in Euclidean metric.

A. Sibson interpolation

Data interpolation based on Voronoï diagrams obtained through different techniques [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF] will describe the scheme proposed by Sibson sake of simplicity, we will consider a cloud known nodal values in a 2-D domain interpolated value at an evaluation point x Initially, the evaluation point x is added to the Voronoï diagram and a new Voronoï cell is created around (Fig. 3). The natural neighbors of x (i.e. the closest nodes) a the nodes that share a Voronoï edge with the Voronoï cell of These will be the nodes used for the interpolation evaluation point. onoï diagrams can be [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF]. In this work, we the scheme proposed by Sibson [START_REF] Sibson | A vector identity for the Dirichlet tessellation[END_REF] influence of each natural neighbor on the evaluation point through Sibson's technique is given by [START_REF] Sibson | A vector identity for the Dirichlet tessellation[END_REF], where each represents the intersection area between the Voronoï cell Voronoï cell related to the , as illustrated by the hashed region in [START_REF] Sibson | A vector identity for the Dirichlet tessellation[END_REF] It is easy to verify that the always positive and less than 1 and that the sum of all te equal to 1. One can equally remark that when n i , the shape function tends towards 1 for node 0 for all other nodes (therefore the

The interpolating character of the Sibson shape functions allows the easy application of the boundary conditions. Furthermore, the shape functions degenerate to a simple linear interpolation at the borders [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] shape function related to a node union of all circumcircles passing through this node.

Last, Fig. 4 Erreur ! Source du renvoi introuvable. shows the elevation of the Sibson shape function on a regular grid of nodes.

B. Galerkin formulation

In order to solve (2), the reduced scalar magnetic potential is used by combine (2), ( 3) and ( 4). This potential is interpolated using [START_REF] Sibson | A vector identity for the Dirichlet tessellation[END_REF]. Then, to the resulted equation and considering permeability in Ω m , we get:

∫ ∫ Γ Ω α α + Ω ϕ α ∇ µ α ∇ m m d b d nj j i mj j i
where φ m is the reduced magnetic potential, induction normal to the boundary and magnetic field normal to the boundary.

In order to simplify the notation, (

[ ] [ ] Q b U S n m = ⋅ + ϕ ⋅ , with , ∫ ∫ Γ Ω α α = Ω α ∇ µ α ∇ = m m U d S j i ij j i ij

IV. BOUNDARY ELEMENT

If no source field in the free space is taken into account, the problem can be represented by a Laplacian equation term of a total scalar magnetic potential solve Laplacian equation transforming the volume integral 2 is easy to verify that the Sibson shape functions (6) are always positive and less than 1 and that the sum of all terms is remark that when x tends towards tends towards 1 for node n i and towards refore the function is interpolant). interpolating character of the Sibson shape functions allows the easy application of the boundary conditions. Furthermore, the shape functions degenerate to a simple linear [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. Fig. 4 (a) depicts the Sibson node. Its support is given by the union of all circumcircles passing through this node.

! Source du renvoi introuvable.(b) the elevation of the Sibson shape function on a regular

∫ Γ Γ µ α = Γ m d H d n i 0 , (7)
is the reduced magnetic potential, b n is the magnetic induction normal to the boundary and H 0n is the source magnetic field normal to the boundary.

In order to simplify the notation, ( 7) can be rewritten as [START_REF] Eibert | Calculation of unbounded field problems in free space by a 3D FEM/BEM-hybrid approach[END_REF] .

0 ∫ Γ Γ µ α = Γ m d H Q and d n i ij j (9)
LEMENT METHOD (BEM) in the free space is taken into account, the problem can be represented by a Laplacian equation in term of a total scalar magnetic potential. The main idea is to equation transforming the volume integral equation over Ω 0 into a surface integral at the boundary Γ m . The classical boundary integral equation is obtained by using the third Green's identity [START_REF] Brebbia | The Boundary Element Method For Engineers[END_REF].

Since we have non-meshed coil in the air, the reduced magnetic potential φ Γ will be introduced. According to (3), we get:

( ) ( ) ( ) ( ) , ) 1 ( , , 0 Γ         ⋅ - - ϕ ∇ = ϕ ∫ Γ Γ Γ d µ G G c m n n h b x x x x
x x 0 0 0 0 0 [START_REF] Ma | A general algorithm for accurate computation of field variables and its derivatives near the boundary in BEM[END_REF] where G(x 0 ,x) is the Green's function ln(1/r), r is the distance between x 0 and a point x on the boundary and c(x 0 ) is the solid angle subtended by the boundary Γ m at x 0 .

The interpolation of the potential and the magnetic induction normal to the boundary are given by:

( ) ( ) j N j j w Γ Γ ϕ = ϕ ∑ Γ x x , ( ) ( ) j n N j j n b w b x x ∑ Γ = , ( 11
)
where N Γ the number of facets of the boundary and w j (x) = 1 if x belongs to the boundary element j and w j (x) = 0 otherwise (0-order). Equation ( 11) can be rewritten, as follow:

[ ] [ ] [ ] n n Γ h T b T H 0 1 = + ϕ . (12)
T, T 1 and H matrices are defined as follows:

( ) Γ = ∫ Γ d G µ T j m ij j i x x , 1 0 , ( ) Γ = ∫ Γ d G T j m ij j i x x , 1 , ( ) Γ ∇ - = ∫ Γ d G c H j m n ii ij j i x x , , (13)
where c ii is equals π, since the collocation point is located on the centroid of each boundary element. There are many ways to perform integration of (13). In this work the well-known analytical expression found in [START_REF] Wang | Analytical evaluation and application of the singularities in boundary element method[END_REF] has been used.

V. COUPLING BOTH METHODS

By considering conservation of the normal component of the magnetic induction and the unicity of the reduced magnetic potential on boundary, we get [START_REF] Eibert | Calculation of unbounded field problems in free space by a 3D FEM/BEM-hybrid approach[END_REF]:

S 0 m ϕ Q U Γ ϕ = . ( 14
) 0 H T n b n 1 h T 0 ⋅
It's important to remark that the BEM portion of the matrix is fully populated and the NEM portion is banded. Moreover the equation ( 14) leads to a linear system which will be solved directly.

By solving the system above, we get the reduced scalar magnetic potential on all nodes and b n on the boundary [START_REF] Meunier | Hybrid finite element boundary element solutions for three dimensional scalar potential problems[END_REF].

VI. RESULTS

The validation of the presented scheme was performed through the simulation of two test cases.

A. Results on a test case with analytical solution

Fig. 5 shows a homogeneous infinite cylinder under the influence of a constant magnetic field (H 0 = 1 A/m). The magnetic domain (cylinder) is modeled using NEM and the unbounded domain (air) using BEM. Fig. 5. Infinite cylinder with a constant permeability immersed in a constant field H0. The total field is evaluated on the path. Fig. 6. Magnitude of the total magnetic field on path (Fig. 5) using the cylinder problem for both analytical and NEM-BEM approach.

The problem was solved with a relatively coarse discretization, generating 232 degrees of freedom (which 80 were on the boundary).Once the problem is solved, the integral equation ( 11) is used to evaluate the reduced scalar magnetic potential at the path showed in Fig. 5. Then, since the source field is known, by ( 3) and (4) the total magnetic field is evaluated. More details on the evaluating of the total magnetic field can be found in [START_REF] Ma | A general algorithm for accurate computation of field variables and its derivatives near the boundary in BEM[END_REF].

Fig. 6 shows the comparison between the analytical result and the NEM-BEM one. A good agreement between the solutions is verified.

B. Results on a real case

Fig. 7 shows a magnetic circuit with a non-meshed coil. In order to test the robustness of the proposed scheme the magnetic induction was computed on two paths: the first one in the gap between the parts of the device and the second one crossing the magnetic circuit.

Once there is no exact analytical solution to this problem, a 2 nd order FEM solution was used as reference. This solution was computed with a very fine mesh. Three different levels of discretization where tested in order to test the convergence of the proposed method. Table 1 shows the number of degrees of freedom (DoF) for each case.

Fig. 8 presents the comparison between the tested cases and the reference. Results show a fast convergence to the reference value as the number of DoF increases. One can remark that the number of DoF required to get an accurate result is much lower than the one of the reference as shown in table I. Fig. 9 illustrates the magnetic induction on the path 2 obtained by NEM-BEM and FEM-BEM. Let us notice that both solutions use the same number of DoF (4,512) and exactly the same node distribution.

A good agreement between the two methods and the reference is achieved. However, one can observe that the NEM-BEM provides a much smoother approximation when compared to FEM-BEM.

VII. CONCLUSION

In this paper, an original hybrid natural element methodboundary element method formulation has been presented. This scheme is very useful, once the flux linkage is automatically taken into account and the mesh dependency is drastically reduced.

The NEM-BEM was able to provide accurate solutions even with very coarse discretizations. Besides that, this new formulation is attractive because of the easiness on the treatment of the unbounded regions and the smoothness of the interpolation inside magnetic regions.

Fig. 1 .

 1 Fig. 1. Magnetostatic problem: bounded and unbounded regions.

Fig. 2

 2 illustrates a Voronoï cell related to a node n i in a 2-D space. The Delaunay triangulation, which is the dual of the Voronoï diagram, is constructed by connecting the nodes whose Voronoï cells have common boundaries. These nodes are called natural neighbors. Alternatively, the vertices of the T Manuscript received June 6, 2014. Corresponding author: D. Martins Araujo (douglas.martins-araujo@fr.schneider-electric.com Voronoï cell are the center of the circumcircles of the Delaunay triangles.

Fig. 2 .

 2 Fig. 2. Voronoï diagram in blue, its associated Delaunay triangulation in red and the circumcircles associated to the Delaunay triangles in grey. Around its Voronoï cell and its 6 natural neighbors, from n1 to

Fig. 3 .

 3 Fig. 3. Sibson natural neighbors and natural coordinates (original Voronoï diagram in blue, natural neighbors in green, additional Voronoï cell due to in black and intersection area between Voronoï cells ofThe influence of each natural neighbor on the evaluation point through Sibson's technique is given by A i (x) represents the intersection area between the Voronoï cell centered on x and the "original" Voronoï cell related to the natural neighbor n i , as illustrated by the hashed region in Fig.3.( ) ( ) ( ) ∑ =

  in 1980. For cloud of nodes with domain, and look for the . added to the "original" Voronoï diagram and a new Voronoï cell is created around it (i.e. the closest nodes) are with the Voronoï cell of x. nodes used for the interpolation at the . Sibson natural neighbors and natural coordinates (original Voronoï ditional Voronoï cell due to x in black and intersection area between Voronoï cells of x and n2 in green).

Fig. 4 .

 4 Fig. 4. (a) Support for NEM shape function and (b) Plot of the Sibson shape function.

Fig. 7 .

 7 Fig. 7. U shaped contactor. Material with linear permeability and a nonmeshed coil. Two different paths have been considered.

Fig. 8 .

 8 Fig.8. Comparison between FEM and the approach presented. Vertical component of the magnetic field on path 1.

Fig. 9 .

 9 Fig.9. Comparison between FEM and the approach presented. Vertical component of the magnetic field on path 2.

TABLE 1

 1

		DEGREES OF FREEDOM
		FEM 2nd	NEM-BEM	NEM-BEM	NEM-BEM
		Reference	Dense	Medium	Coarse
	DOF	129,751	3,745	265	77

. (a) Support for NEM shape function and (b) Plot of the Sibson shape the reduced scalar magnetic potential combine (2), (3) and (4). This potential is by applying the Galerkin method the resulted equation and considering a constant