2.1. Admissible trajectories modeling

We start with definitions.

When optimizing a trajectory with respect to a given criterion, the initial and final states are often constrained, that is to say the optimization is performed in an affine subspace modeling these endpoints conditions. This subspace is now introduced.

Definition 2 (Endpoints conditions). Let $ {y}_0,{y}_T\in {\unicode{x211D}}^D $ . We define the set $ \mathcal{D}\left({y}_0,{y}_T\right)\subset \mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right) $ as

$$ y\in \mathcal{D}\left({y}_0,{y}_T\right)\hskip2em \iff \hskip2em \left\{\begin{array}{l}y(0)={y}_0,\\ {}y(T)={y}_T.\end{array}\right. $$

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 user-defined constraints. We define now the set of trajectories verifying such additional constraints.

Definition 3 (Additional constraints). For $ \mathrm{\ell}=1,\dots, L $ , let $ {g}_{\mathrm{\ell}} $ be a real-valued function defined on $ {\unicode{x211D}}^D $ . We define the set $ \mathcal{G}\subset \mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right) $ as the set of trajectories over $ \left[0,T\right] $ satisfying the following $ L $ inequality constraints given by the functions $ {g}_{\mathrm{\ell}} $ , that is

$$ y\in \mathcal{G}\hskip2em \iff \hskip2em \forall \mathrm{\ell}=1,\dots, L\hskip1em \forall t\in \left[0,T\right]\hskip2em {g}_{\mathrm{\ell}}\left(y(t)\right)\leqslant 0. $$

Lastly, we introduce the set of admissible trajectories which satisfy both the endpoints conditions and the additional constraints.

Definition 4 (Admissible trajectory). We define the set $ {\mathcal{A}}_{\mathcal{G}}\left({y}_0,{y}_T\right)\subset \mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right) $ as follows:

$$ {\mathcal{A}}_{\mathcal{G}}\left({y}_0,{y}_T\right):= \mathcal{D}\left({y}_0,{y}_T\right)\cap \mathcal{G}. $$

Any element of $ {\mathcal{A}}_{\mathcal{G}}\left({y}_0,{y}_T\right) $ will be called an admissible trajectory.

2.2. Projection for a finite-dimensional optimization problem

In our approach, a theoretical optimization 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 (i.e., the set of projected trajectories in Figure 1).

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 5. Let $ {\left\{{\varphi}_k\right\}}_{k=1}^{+\infty } $ be an orthonormal basis of $ {L}^2\left(\left[0,T\right],\unicode{x211D}\right) $ ^{Footnote 2} with respect to the inner product

$$ \left\langle f,g\right\rangle ={\int}_0^Tf(t)\hskip0.1em g(t)\hskip0.1em dt, $$

such that each $ {\varphi}_k $ is continuous on $ \left[0,T\right] $ and let $ \mathcal{K}:= {\left\{{K}_d\right\}}_{d=1}^D $ be a sequence of integers with $ K:= {\sum}_{d=1}^D{K}_d $ . We define the space of projected trajectories $ {\mathcal{Y}}_{\mathcal{K}}\left(0,T\right)\subset \mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right) $ over $ \left[0,T\right] $ as

$$ {\mathcal{Y}}_{\mathcal{K}}\left(0,T\right):= \prod \limits_{d=1}^D span{\left\{{\varphi}_k\right\}}_{k=1}^{K_d}. $$

If there is no risk of confusion, we write $ {\mathcal{Y}}_{\mathcal{K}}:= {\mathcal{Y}}_{\mathcal{K}}\left(0,T\right) $ for the sake of readability.

Remark 1. From the above definition, any projected trajectory $ y\in {\mathcal{Y}}_{\mathcal{K}} $ is associated with a unique vector

$$ c={\left({c}_1^{(1)},\dots, {c}_{K_1}^{(1)},{c}_1^{(2)},\dots, {c}_{K_2}^{(2)},\dots, {c}_1^{(D)},\dots, {c}_{K_D}^{(D)}\right)}^T\in {\unicode{x211D}}^K $$

defined by

(2) $$ {c}_k^{(d)}:= \left\langle {y}^{(d)},{\varphi}_k\right\rangle ={\int}_0^T{y}^{(d)}(t)\hskip0.1em {\varphi}_k(t)\hskip0.1em dt. $$

In other words, the vector $ c $ is the image of the trajectory $ y $ by the projection operator $ \Phi :\mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right)\to {\unicode{x211D}}^K $ defined by $ \Phi y:= c $ , whose restriction $ {\left.\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}} $ is bijective (as the Cartesian product of bijective operators). In particular, the spaces $ {\mathcal{Y}}_{\mathcal{K}} $ and $ {\unicode{x211D}}^K $ are isomorphic, that is $ {\mathcal{Y}}_{\mathcal{K}}\simeq {\unicode{x211D}}^K $ .

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

Proposition 1. A trajectory $ y\in {\mathcal{Y}}_{\mathcal{K}} $ belongs to $ \mathcal{D}\left({y}_0,{y}_T\right) $ if and only if its associated vector $ c:= \Phi y\in {\unicode{x211D}}^K $ satisfies the linear system

(3) $$ A\left(0,T\right)\hskip0.1em c=\Gamma, $$

where the matrix $ A\left(0,T\right)\in {\unicode{x211D}}^{2D\times K} $ and the vector $ \Gamma \in {\unicode{x211D}}^{2D} $ are given by

$$ A\left(0,T\right):= \left(\begin{array}{ccccccc}{\varphi}_1(0)& \dots & {\varphi}_{K_1}(0)& & & & \\ {}& & & \ddots & & & \\ {}& & & & {\varphi}_1(0)& \dots & {\varphi}_{K_D}(0)\\ {}{\varphi}_1(T)& \dots & {\varphi}_{K_1}(T)& & & & \\ {}& & & \ddots & & & \\ {}& & & & {\varphi}_1(T)& \dots & {\varphi}_{K_D}(T)\end{array}\right)\hskip1em ,\hskip1em \Gamma := \left(\begin{array}{c}{y}_0\\ {}{y}_T\end{array}\right). $$

Proof. Let $ y\in {\mathcal{Y}}_{\mathcal{K}} $ and let $ c:= \Phi y\in {\unicode{x211D}}^K $ . By the definition of the matrix $ A\left(0,T\right) $ , we have

$$ {\displaystyle \begin{array}{c}A\left(0,T\right)\hskip0.1em c=A\left(0,T\right){\left({c}_1^{(1)},\dots, {c}_{K_1}^{(1)},{c}_1^{(2)},\dots, {c}_{K_2}^{(2)},\dots, {c}_1^{(D)},\dots, {c}_{K_D}^{(D)}\right)}^T\\ {}\hskip18.9em ={\left(\sum \limits_{k=1}^{K_1}{c}_k^{(1)}{\varphi}_k(0),\dots, \sum \limits_{k=1}^{K_D}{c}_k^{(D)}{\varphi}_k(0),\dots, \sum \limits_{k=1}^{K_1}{c}_k^{(1)}{\varphi}_k(T),\dots, \sum \limits_{k=1}^{K_D}{c}_k^{(D)}{\varphi}_k(T)\right)}^T\\ {}\hskip-21.5em =\left(\begin{array}{c}y(0)\\ {}y(T)\end{array}\right).\end{array}} $$

The conclusion follows directly from the preceding relation.

2.3. Reference trajectories modeling

Let us now suppose that we have access to $ I $ recorded trajectories $ {y}_{R_1},\dots, {y}_{R_I} $ , called reference trajectories, coming from some experiments. We propose here an example of a statistical modeling for these reference trajectories, permitting especially to exhibit some linear properties. This modeling will allow to take advantage of the information contained in these recorded trajectories when deriving optimization problems.

These trajectories being recorded, they are in particular admissible and we assume that they belong to the space $ {\mathcal{Y}}_{\mathcal{K}}\left(0,T\right) $ . 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}\in {\unicode{x211D}}^K $ . Moreover, we consider each reference trajectory as a noisy observation of a certain admissible and projected trajectory $ {y}_{\ast } $ . In other words, we suppose that there exists a trajectory $ {y}_{\ast}\in {\mathcal{Y}}_{\mathcal{K}}\cap {\mathcal{A}}_{\mathcal{G}}\left({y}_0,{y}_T\right) $ associated with a vector $ {c}_{\ast}\in {\unicode{x211D}}^K $ satisfying

$$ \forall i=1,\dots, I\hskip2em {c}_{R_i}={c}_{\ast }+{\varepsilon}_i. $$

The noise $ {\varepsilon}_i $ is here assumed to be a centered GaussianFootnote ³ whose covariance matrix $ {\Sigma}_i $ is of the form

$$ {\Sigma}_i=\frac{1}{2{\omega}_i}\hskip0.1em \Sigma, $$

where $ \Sigma \in {\unicode{x211D}}^{K\times 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 $ {\omega}_i $ are here considered as weights so we require $ {\sum}_{i=1}^I{\omega}_i=1 $ ; each $ {\omega}_i $ plays actually the role of a noise intensity. Further from the hypothesis that the trajectory $ {y}_{\ast } $ and all the reference trajectories $ {y}_{R_i} $ verify the same endpoints conditions, we deduce

$$ A\hskip0.1em {c}_{R_i}=A\hskip0.1em {c}_{\ast }+A\hskip0.1em {\varepsilon}_i\hskip2em \iff \hskip2em A\hskip0.1em {\varepsilon}_i={0}_{{\unicode{x211D}}^{2D}}\hskip2em \iff \hskip2em {\varepsilon}_i\in \ker A, $$

for all $ i=1,\dots, I $ (we shorten $ A\left(0,T\right) $ in $ A $ when the context is clear). Hence, the reference vector $ {c}_{\ast } $ satisfies the following $ I $ systems:

(4) $$ \left\{\begin{array}{l}{c}_{R_i}={c}_{\ast }+{\varepsilon}_i,\\ {}{\varepsilon}_i\sim \mathcal{N}\left({0}_{{\unicode{x211D}}^K},{\Sigma}_i\right)\\ {}{\varepsilon}_i\in \ker A.\end{array}\right., $$

To establish a more explicit system which is equivalent to the preceding one, we require the following preliminary proposition. Here, we diagonalize the matrices $ \Sigma $ and $ {A}^TA $ 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 modeling: the endpoints conditions modelled by $ A $ imply linear relations within the components of the vectors, which should be reflected by the covariance matrix $ \Sigma $ . The following result will be helpful to establish the upcoming proposition 3.

Proposition 2. We define $ \sigma := \mathit{\operatorname{rank}}\hskip0.1em \Sigma $ and $ a:= \mathit{\operatorname{rank}}\hskip0.1em {A}^TA $ . In the setting of system (4), we have $ \sigma +a\leqslant K $ and there exist an orthogonal matrix $ V\in {\unicode{x211D}}^{K\times K} $ and two matrices $ {\Lambda}_{\Sigma}\in {\unicode{x211D}}^{K\times K} $ and $ {\Lambda}_A\in {\unicode{x211D}}^{K\times K} $ of the following form:

$$ {\Lambda}_{\Sigma}=\left(\begin{array}{cc}{\Lambda}_{\Sigma, 1}& {0}_{{\unicode{x211D}}^{\sigma \times \left(K-\sigma \right)}}\\ {}{0}_{{\unicode{x211D}}^{\left(K-\sigma \right)\times \sigma }}& {0}_{{\unicode{x211D}}^{\left(K-\sigma \right)\times \left(K-\sigma \right)}}\end{array}\right)\hskip2em ,\hskip2em {\Lambda}_A=\left(\begin{array}{cc}{0}_{{\unicode{x211D}}^{\left(K-a\right)\times \left(K-a\right)}}& {0}_{{\unicode{x211D}}^{\left(K-a\right)\times a}}\\ {}{0}_{{\unicode{x211D}}^{a\times \left(K-a\right)}}& {\Lambda}_{A,2}\end{array}\right), $$

where $ {\Lambda}_{\Sigma, 1}\in {\unicode{x211D}}^{\sigma \times \sigma } $ and $ {\Lambda}_{A,2}\in {\unicode{x211D}}^{a\times a} $ are diagonal matrices with positive elements, such that

$$ \Sigma =V{\Lambda}_{\Sigma}{V}^T\hskip2em ,\hskip2em {A}^TA=V{\Lambda}_A{V}^T. $$

Proof. The proof starts by noticing

(5) $$ \Sigma \hskip0.1em {A}^TA={A}^TA\hskip0.1em \Sigma ={0}_{{\unicode{x211D}}^{K\times K}}. $$

Indeed using the hypothesis $ {\varepsilon}_i\in \ker \hskip0.5em A $ for any $ i=1,\dots, I $ gives

$$ \Sigma \hskip0.1em {A}^TA=2{\omega}_i\hskip0.1em {\Sigma}_i\hskip0.1em {A}^TA=2{\omega}_i\hskip0.1em \unicode{x1D53C}\left({\varepsilon}_i{\varepsilon}_i^T\right)\hskip0.1em {A}^TA=2{\omega}_i\hskip0.1em \unicode{x1D53C}\left({\varepsilon}_i\hskip0.1em {\left(A{\varepsilon}_i\right)}^T\right)\hskip0.1em A={0}_{{\unicode{x211D}}^{K\times K}}; $$

similar arguments prove the second equality in (5). First, we can deduce

(6) $$ \operatorname{Im}\hskip0.3em \Sigma \hskip0.3em \subseteq \hskip0.3em \ker \hskip0.3em {A}^TA, $$

which leads to $ \sigma \leqslant K-a $ by the rank-nullity theorem. Equalities (5) show also that $ \Sigma $ and $ {A}^TA $ are simultaneously diagonalizable (since they commute) so there exists an orthogonal matrix $ V\in {\unicode{x211D}}^{K\times K} $ such that

(7) $$ \Sigma =V{\Lambda}_{\Sigma}{V}^T\hskip2em ,\hskip2em {A}^TA=V{\Lambda}_A{V}^T, $$

where $ {\Lambda}_{\Sigma}\in {\unicode{x211D}}^{K\times K} $ and $ {\Lambda}_A\in {\unicode{x211D}}^{K\times K} $ are diagonal matrices. Permuting if necessary columns of $ V $ , we can write the matrix $ {\Lambda}_{\Sigma} $ as follows:

(8) $$ {\Lambda}_{\Sigma}=\left(\begin{array}{cc}{\Lambda}_{\Sigma, 1}& {0}_{{\unicode{x211D}}^{\sigma \times \left(K-\sigma \right)}}\\ {}{0}_{{\unicode{x211D}}^{\left(K-\sigma \right)\times \sigma }}& {0}_{{\unicode{x211D}}^{\left(K-\sigma \right)\times \left(K-\sigma \right)}}\end{array}\right); $$

in other words, the $ \sigma $ first column vectors of $ V $ span the image of $ \Sigma $ . From the inclusion (6), we deduce that these vectors belong to the null space of $ {A}^TA $ . Hence, the $ \sigma $ first diagonal elements of $ {\Lambda}_A $ are equal to zero and, up to a permutation of the $ K-\sigma $ last column vectors of $ V $ , we can write

$$ {\Lambda}_A=\left(\begin{array}{cc}{0}_{{\unicode{x211D}}^{\left(K-a\right)\times \left(K-a\right)}}& {0}_{{\unicode{x211D}}^{\left(K-a\right)\times a}}\\ {}{0}_{{\unicode{x211D}}^{a\times \left(K-a\right)}}& {\Lambda}_{A,2}\end{array}\right), $$

which ends the proof. □

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

$$ \operatorname{Im}\hskip0.5em {A}^TA\subseteq \ker \hskip0.3em \Sigma, $$

showing that $ \Sigma $ is singular. Consequently the Gaussian noise $ {\varepsilon}_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 using the preceding result which allows to decompose the space $ {\unicode{x211D}}^K $ into three orthogonal subspaces. We prove that the restriction of the noise $ {\varepsilon}_i $ to the first subspace is a non-degenerate 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 $ \Sigma $ .

Prior to this, let us write the matrix $ V\in {\unicode{x211D}}^{K\times K} $ introduced in Proposition 2 as follows:

$$ V=\left({V}_1\hskip1em {V}_2\hskip1em {V}_3\right), $$

where $ {V}_1\in {\unicode{x211D}}^{K\times \sigma } $ , $ {V}_2\in {\unicode{x211D}}^{K\times K-\sigma -a} $ and $ {V}_3\in {\unicode{x211D}}^{K\times a} $ . We emphasize that the column-vectors of the matrices $ {V}_1 $ and $ {V}_3 $ do not overlap according to the property $ \sigma +a\leqslant K $ proved in proposition 2. In particular, the matrix $ {V}_2 $ has to be considered only in the case $ \sigma +a<K $ . Further for any $ c\in {\unicode{x211D}}^K $ , we will use the notation

$$ \tilde{c}:= {V}^Tc\hskip2em ,\hskip2em {\tilde{c}}_{\mathrm{\ell}}:= {V}_{\mathrm{\ell}}^Tc, $$

for $ \mathrm{\ell}=\mathrm{1,2,3} $ . Finally, we consider the singular value decomposition of $ A $ coming from the diagonalization of the symmetric matrix $ {A}^TA $ with $ V $ :

$$ A={US}_A{V}^T, $$

where $ U\in {\unicode{x211D}}^{2D\times 2D} $ is orthogonal and $ {S}_A\in {\unicode{x211D}}^{2D\times K} $ is a rectangular diagonal matrix of the following form:

(9) $$ {S}_A=\left({0}_{{\unicode{x211D}}^{2D\times K-2D}}\hskip1em {S}_{A,2}\right), $$

with $ {S}_{A,2}:= \sqrt{\Lambda_{A,2}}\in {\unicode{x211D}}^{2D\times 2D} $ .

Proposition 3. Suppose that the matrix $ A $ is full rank, that is $ a=2D $ . Then for any $ i=1,\dots, I $ , system (4) is equivalent to the following one:

(10) $$ \left\{\begin{array}{l}{\tilde{c}}_{R_i,1}={\tilde{c}}_{\ast, 1}+{\tilde{\varepsilon}}_{i,1},\\ {}{\tilde{\varepsilon}}_{i,1}\sim \mathcal{N}\left({0}_{{\unicode{x211D}}^{\sigma }},\frac{1}{2{\omega}_i}\hskip0.1em {\Lambda}_{\Sigma, 1}\right)\\ {}{\tilde{c}}_{\ast, 2}={V}_2^T{c}_{R_i},\\ {}{\tilde{c}}_{\ast, 3}={S}_{A,2}^{-1}\hskip0.1em {U}^T\Gamma .\end{array},\right. $$

Proof. We first prove that system (4) is equivalent to

(11) $$ \left\{\begin{array}{l}{\tilde{c}}_{R_i}={\tilde{c}}_{\ast }+{\tilde{\varepsilon}}_i,\\ {}{\tilde{\varepsilon}}_i\sim \mathcal{N}\left({0}_{{\unicode{x211D}}^K},\frac{1}{2{\omega}_i}\hskip0.1em {\Lambda}_{\Sigma}\right)\\ {}{S}_A{\tilde{c}}_{\ast }={U}^T\Gamma .\end{array}\right., $$

The matrix $ V $ being orthogonal, it is nonsingular and so we have for all $ i=1,\dots, I $ ,

$$ {c}_{R_i}={c}_{\ast }+{\varepsilon}_i\hskip2em \iff \hskip2em {\tilde{c}}_{R_i}={\tilde{c}}_{\ast }+{\tilde{\varepsilon}}_i, $$

and, since $ {\Sigma}_i=\frac{1}{2{\omega}_i}\hskip0.1em \Sigma =\frac{1}{2{\omega}_i}\hskip0.1em V{\Lambda}_{\Sigma}{V}^T $ , we obtain

$$ {\varepsilon}_i\sim \mathcal{N}\left({0}_{{\unicode{x211D}}^K},{\Sigma}_i\right)\hskip2em \iff \hskip2em {\tilde{\varepsilon}}_i\sim \mathcal{N}\left({0}_{{\unicode{x211D}}^K},\frac{1}{2{\omega}_i}\hskip0.1em {\Lambda}_{\Sigma}\right). $$

Finally, the property $ {\varepsilon}_i\in \ker \hskip0.3em A $ is equivalent to

$$ A\hskip0.1em {c}_{\ast }=\Gamma \hskip2em \iff \hskip2em {US}_A{V}^T{c}_{\ast }=\Gamma \hskip2em \iff \hskip2em {S}_A\hskip0.1em {\tilde{c}}_{\ast }={U}^T\Gamma, $$

proving that the systems (4) and (11) are equivalent. Now the fact that the $ K-\sigma $ last diagonal elements of $ {\Lambda}_{\Sigma} $ are zero implies that the components $ {\tilde{c}}_{\ast, 2}\in {\unicode{x211D}}^{K-\sigma -2D} $ and $ {\tilde{c}}_{\ast, 3}\in {\unicode{x211D}}^{2D} $ are constant. From the first equality of (11), we have on one side

$$ {\tilde{c}}_{R_i,2}={\tilde{c}}_{\ast, 2}\hskip2em \iff \hskip2em {V}_2^T{c}_{R_i}={\tilde{c}}_{\ast, 2}, $$

for any $ i=1,\dots, I $ . On the other side, combining the last relation of the system (11) with the form of the matrix $ {S}_A $ given in (9) yields

$$ {\displaystyle \begin{array}{l}{S}_A\hskip0.1em {\tilde{c}}_{\ast }={U}^T\Gamma \hskip2em \iff \hskip2em {S}_{A,2}\hskip0.1em {\tilde{c}}_{\ast, 3}={U}^T\Gamma \\ {}\hskip10em \iff \hskip2em {\tilde{c}}_{\ast, 3}={S}_{A,2}^{-1}\hskip0.1em {U}^T\Gamma, \end{array}} $$

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 nonsingular).□

The above decomposition gives us access to nondegenerated density of $ {\tilde{c}}_{R_i,1} $ given $ {\tilde{c}}_{\ast, 1} $ which is later denoted by $ u\left({\tilde{c}}_{R_i,1}|{\tilde{c}}_{\ast, 1}\right) $ . In next section, we will assume a prior distribution on $ {\tilde{c}}_{\ast, 1} $ with high density for low values of the cost function $ F $ .

2.4. A trajectory optimization problem via a MAP approach

Before introducing the Bayesian framework, let first recall that we are interested in minimizing a certain cost function $ F:\mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right)\to \unicode{x211D} $ over the set of projected and admissible trajectories $ {\mathcal{Y}}_{\mathcal{K}}\cap {\mathcal{A}}_{\mathcal{G}}\left({y}_0,{y}_T\right) $ . As explained previously, we propose here a methodology leading to a constrained optimization problem based on the reference trajectories and designed to provide realistic trajectories (we refer again to Figure 1). 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 modeling of the reference trajectories given in Proposition 3 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 $ {\mathcal{Y}}_{\mathcal{K}} $ which is isomorphic to $ {\unicode{x211D}}^K $ . So each trajectory is here described by its associated vector in $ {\unicode{x211D}}^K $ , permitting in particular to define distributions over finite-dimensional spaces. We also recall that the reference trajectories are interpreted as noisy observations of a certain $ {y}_{\ast } $ associated with a $ {c}_{\ast } $ . According to Proposition 3, this vector complies with some affine conditions which are described by the following subspace $ {\mathcal{V}}_1 $ :

(12) $$ c\in {\mathcal{V}}_1\hskip2em \iff \hskip2em \left\{\begin{array}{l}{V}_2^Tc={V}_2^T{c}_{R_i},\\ {}{V}_3^Tc={S}_{A,2}^{-1}\hskip0.1em {U}^T\Gamma .\end{array}\right. $$

Hence, a vector $ c $ belonging to $ {\mathcal{V}}_1 $ is described only through its component $ {\tilde{c}}_1:= {V}_1^Tc $ . In addition, we note that the definition of $ {\mathcal{V}}_1 $ does not depend actually on the choice of $ i $ since $ {V}_2^T{c}_{R_i} $ has been proved to be constant in Proposition 3. Further, we emphasize that the matrix $ A $ is supposed to be full rank in this case and we have $ {\mathcal{V}}_1\simeq {\unicode{x211D}}^{\sigma } $ ; we recall that $ \sigma $ is the rank of the covariance matrix $ \Sigma $ .

Let us now define the cost function $ F $ over the spaces $ {\unicode{x211D}}^K $ and $ {\mathcal{V}}_1 $ . This is necessary to define the prior distribution and to establish our optimization problem.

Remark 3. From the previous definition, we observe that for any $ y\in {\mathcal{Y}}_{\mathcal{K}} $ and its associated vector $ c\in {\unicode{x211D}}^K $ , we have

$$ {\left.\overset{\check{}}{F}(c)=F\Big(\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\Big)=F(y). $$

Further for any $ c\in {\mathcal{V}}_1 $ , we have

$$ \overset{\check{}}{F}(c)=F\left(\Phi {\left|{}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\left)=F\right(\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}V\tilde{c}\right)=F\left({\left.\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}\hskip0.1em V{\left({\tilde{c}}_1^T\hskip1em {c}_{R_i}^T{V}_2\hskip1em {\Gamma}^TU\hskip0.1em {\left({S}_{A,2}^{-1}\right)}^T\right)}^T\right)=\tilde{F}\left({\tilde{c}}_1\right). $$

We deduce that $ \tilde{F} $ is actually the restriction of $ \overset{\check{}}{F} $ to the subspace $ {\mathcal{V}}_1 $ .

From now on, the trajectory $ {y}_{\ast } $ and the associated vector $ {c}_{\ast } $ will be considered as random variables and will be denoted by $ y $ and $ c $ . We are interested in the posterior distribution

$$ u\left({\tilde{c}}_1\hskip0.1em |\hskip0.1em {\tilde{c}}_{R_1,1},\dots, {\tilde{c}}_{R_I,1}\right), $$

which depends only on the free component $ {\tilde{c}}_1 $ of $ c\in {\mathcal{V}}_1 $ , the two other ones $ {\tilde{c}}_2 $ and $ {\tilde{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\left({\tilde{c}}_1\hskip0.1em |\hskip0.1em {\tilde{c}}_{R_1,1},\dots, {\tilde{c}}_{R_I,1}\right)\propto u\left({\tilde{c}}_{R_1,1},\dots, {\tilde{c}}_{R_I,1}\hskip0.1em |\hskip0.1em {\tilde{c}}_1\right)\hskip0.1em u\left({\tilde{c}}_1\right). $$

Assuming now that the vectors $ {\tilde{c}}_{R_i,1} $ are independent gives

$$ u\left({\tilde{c}}_{R_1,1},\dots, {\tilde{c}}_{R_I,1}\hskip0.1em |\hskip0.1em {\tilde{c}}_1\right)\hskip0.1em u\left({\tilde{c}}_1\right)=\prod \limits_{i=1}^Iu\left({\tilde{c}}_{R_i,1}\hskip0.1em |\hskip0.1em {\tilde{c}}_1\right)\hskip0.1em u\left({\tilde{c}}_1\right). $$

The above likelihood is given by the modeling of the reference trajectories detailed in Proposition 3. In this case, we have

$$ u\left({\tilde{c}}_{R_i,1}\hskip0.1em |\hskip0.1em {\tilde{c}}_1\right)\propto \exp \left(-{\omega}_i{\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)}^T{\Lambda}_{\Sigma, 1}^{-1}\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)\right). $$

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:

(13) $$ u\left({\tilde{c}}_1\right)\propto \exp \left(-{\kappa}^{-1}\tilde{F}\left({\tilde{c}}_1\right)\right), $$

where $ \kappa >0 $ . Putting everything together and taking the negative of the logarithm gives the following minimization problem, whose solution is the MAP estimator:

(14) $$ \left\{\begin{array}{l}{\tilde{c}}_1^{\star}\in \arg \hskip0.1em \underset{{\tilde{c}}_1\in {\unicode{x211D}}^{\sigma }}{\min}\tilde{F}\left({\tilde{c}}_1\right)+\kappa \sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)}^T{\Lambda}_{\Sigma, 1}^{-1}\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right),\\ {}{\tilde{c}}_2={V}_2^T{c}_{R_i},\\ {}{\tilde{c}}_3={S}_{A,2}^{-1}\hskip0.1em {U}^T\Gamma .\end{array}\;\right. $$

where $ i $ is arbitrarily chosen in $ \left\{1,\dots, I\right\} $ .

Let us now rewrite the above optimization problem with respect to the variable $ c=V\tilde{c}\in {\unicode{x211D}}^K $ in order to make it more interpretable.

Proposition 4. The optimization problem (14) is equivalent to the following one:

(15) $$ {c}^{\star}\in \underset{c\in {\mathcal{V}}_1}{\arg \hskip0.1em \min}\;\overset{\check{}}{F}(c)+\kappa \sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left(c-{c}_{R_i}\right)}^T{\Sigma}^{\dagger}\left(c-{c}_{R_i}\right), $$

where $ {\Sigma}^{\dagger}\in {\unicode{x211D}}^{K\times K} $ denotes the pseudoinverse of the matrix $ \Sigma $ .

Proof. From (8), we deduce

$$ {\displaystyle \begin{array}{c}\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)}^T{\Lambda}_{\Sigma, 1}^{-1}\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)=\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left(\tilde{c}-{\tilde{c}}_{R_i}\right)}^T{\Lambda}_{\Sigma}^{\dagger}\left(\tilde{c}-{\tilde{c}}_{R_i}\right)\\ {}\hskip25em =\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left(c-{c}_{R_i}\right)}^T\hskip0.1em V{\Lambda}_{\Sigma}^{\dagger }{V}^T\hskip0.1em \left(c-{c}_{R_i}\right)\\ {}\hskip21.5em =\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left(c-{c}_{R_i}\right)}^T{\Sigma}^{\dagger}\left(c-{c}_{R_i}\right).\end{array}} $$

And from the proof of proposition 3, we have

$$ A\hskip0.1em c=\Gamma \hskip2em \iff \hskip2em {\tilde{c}}_3={S}_{A,2}^{-1}\hskip0.1em {U}^T\Gamma, $$

proving that $ c\in {\mathcal{V}}_1 $ .

To conclude, let us comment on this optimization problem.

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

   (16) $$ {\displaystyle \begin{array}{ll}& \underset{{\tilde{c}}_1\in {\unicode{x211D}}^{\sigma }}{\min}\tilde{F}\left({\tilde{c}}_1\right)\\ {}\mathrm{s}.\mathrm{t}.& \sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)}^T{\Lambda}_{\Sigma, 1}^{-1}\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)\leqslant \tilde{\kappa}\end{array}}, $$

   where $ \lambda \geqslant 0 $ . Here, we suppose that $ \tilde{F} $ is strictly convex and that the problem (16) has a solution (which is then unique). By Slater’s theorem (Boyd and Vandenberghe, Reference Boyd and Vandenberghe2004, Section 5.2.3), the strong duality holds for the problem (16). It can then be proved that there exists a certain $ {\lambda}^{\star}\geqslant 0 $ such that the solution of (16) is the minimizer of the strictly convex function

   $$ {\tilde{c}}_1\mapsto \tilde{F}\left({\tilde{c}}_1\right)+{\lambda}^{\star}\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right)}^T{\Lambda}_{\Sigma, 1}^{-1}\left({\tilde{c}}_1-{\tilde{c}}_{R_i,1}\right), $$

   which is actually the objective function of the optimization problem (14) for $ \kappa ={\lambda}^{\star } $ . Hence, the problem (14) minimizes the cost $ \tilde{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 $ \tilde{F} $ then one could expect the solution of (14) to be equal to this optimal trajectory.

2. 2. Further the optimization problem (15) takes into account the endpoints conditions through the subspace $ {\mathcal{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 $ \kappa >0 $ . We refer to Section 2.6 for more details on an iterative method for the tuning of $ \kappa $ . In particular, a right choice for this parameter is expected to provide an optimized trajectory with a realistic behavior. This is for instance illustrated in Section 4.

3. 3. Taking into account the linear information from the available data through the covariance matrix $ \Sigma $ allows to restrict the search to the subspace $ {\mathcal{V}}_1 $ describing these relations. This is of particular interest when implicit relations (modeled by the submatrix $ {V}_2 $ ) are revealed by the estimation of $ \Sigma $ on the basis of the reference trajectories; in this case, these implicit relations may not be known by the expert.

4. 4. The optimization problem (15) has linear constraints and a quadratic penalized term. For instance, if the cost function $ \overset{\check{}}{F} $ is a convex function then we obtain a convex problem for which efficient algorithms exist.

2.5. Quadratic cost for a convex optimization 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, that is

(17) $$ \forall y\in \mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^D\right)\hskip2em F(y)={\int}_0^Tf\left(y(t)\right)\hskip0.1em dt, $$

where $ f:{\unicode{x211D}}^D\to \unicode{x211D} $ is quadratic. Even though such a setting may appear to be restrictive, we emphasize that quadratic models may lead to highly accurate approximations of variables, as it is illustrated in Section 4. For a quadratic instantaneous cost, the associated function $ \overset{\check{}}{F}:{\unicode{x211D}}^K\to \unicode{x211D} $ can be proved to be quadratic as well and can be explicitly computed. In the following result, we provide a quadratic optimization problem equivalent to (15).

Proposition 5. Suppose that the cost function $ F $ is of the form (17) with $ f $ quadratic. Then the optimization problem (15) is equivalent to the following one:

(18) $$ {c}^{\star}\in \underset{c\in {\mathcal{V}}_1}{\arg \hskip0.1em \min}\;{c}^T\left(\overset{\check{}}{Q}+\kappa \hskip0.1em {\Sigma}^{\dagger}\right)c+{\left(\overset{\check{}}{w}-2\kappa \sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\Sigma}^{\dagger }{c}_{R_i}\right)}^Tc, $$

where $ \overset{\check{}}{Q}\in {\unicode{x211D}}^{K\times K} $ and $ \overset{\check{}}{w}\in {\unicode{x211D}}^K $ can be explicitly computed from $ f $ .

Proof. We defer the proof to Appendix A.

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

2.6. Iterative process to comply with additional constraints

As explained in Section 2.4, the trajectory optimization 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 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 $ {\left.{y}^{\star }=\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}{c}^{\star } $ , where $ {c}^{\star}\in {\mathcal{V}}_1 $ is the solution of the optimization problem (15), verifies the additional constraints, that is belongs to the set $ \mathcal{G} $ . A first solution would be to add the constraint $ {\left.\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\in \mathcal{G} $ in the optimization problem (15). However, depending on the nature of the constraints functions $ {g}_{\mathrm{\ell}} $ , 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 $ \kappa >0 $ appearing in the problem (15).

First of all, let us factorize the problem (15) by $ \kappa $ to obtain the following new one for the sake of presentation:

(19) $$ {c}^{\star}\in \underset{c\in {\mathcal{V}}_1}{\arg \hskip0.1em \min}\;\nu \hskip0.1em \overset{\check{}}{F}(c)+\sum \limits_{i=1}^I{\omega}_i\hskip0.1em {\left(c-{c}_{R_i}\right)}^T{\Sigma}^{\dagger}\left(c-{c}_{R_i}\right), $$

where $ \nu := {\kappa}^{-1} $ . On one hand, we observe that the solution of the optimization problem (19) for the limit case $ \nu =0 $ is given by $ {\sum}_{i=1}^I{\omega}_i\hskip0.1em {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 optimize the cost function $ F $ . On the other hand, for very large $ \nu >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 $ {\nu}^{\star }>0 $ in order to reach a trade-off between optimizing and remaining close to the reference trajectories to comply with the additional constraints. Many methods can be developed to find such a $ {\nu}^{\star } $ but here we propose a very simple one based on a binary search algorithm. We exploit the assumption that the solution for $ \nu =0 $ (i.e., the reference trajectories average) is admissible. We proceed then as follows:

1. 1. We set firstly a maximal value $ {\nu}_{max}>0 $ so that the solution of (19) with $ {\nu}_{max} $ does not satisfy the constraints. For instance, we choose a large value, solve the optimization problem and check the constraints. If they are still satisfied, we choose a larger $ {\nu}_{max} $ .

2. 2. Then apply a binary search between $ 0 $ and $ {\nu}_{max} $ : solve the optimization problem for $ \frac{\nu_{max}}{2} $ ; if the resulting solution verifies the constraints, solve the problem again for $ \frac{3{\nu}_{max}}{4} $ ; else solve it for $ \frac{\nu_{max}}{4} $ .

3. 3. Continue the process until a user-defined stopping criterion.

Other iterative processes can also be considered.

2.7. Confidence bounds on the integrated cost

In practice the cost function $ F $ considered is an estimation of the true cost $ {F}^{\star } $ , a random variable which cannot be fully predicted based on $ y $ . If the distribution $ F(y) $ was known we could deduce a confidence bound on $ {F}^{\star } $ . This is for instance possible by considering multivariate functional regression Ramsay et al. (Reference Ramsay, Hooker, Campbell and Cao2007).

The simplest case from the estimation point of view is to consider that $ {F}^{\star } $ is the integral of some instantaneous consumption function $ {f}^{\star } $ as in Section 2.5, and to estimate the parameters of the standard multivariate regression

$$ {f}^{\star}\left(y(t)\right)=f\left(y(t)\right)+\varepsilon (t), $$

where the random noise $ \varepsilon (t) $ is assumed to follow a centered Gaussian distribution with variance $ \sigma $ . In this case, $ {F}^{\star } $ can be expressed as the integral of a stochastic process

$$ {F}^{\star }(y):= {\int}_0^T{f}^{\star}\left(y(t)\right)\hskip0.1em dt=F(y)+{\int}_0^T\varepsilon (t)\hskip0.1em dt, $$

then assuming that $ {\left(\varepsilon (t)\right)}_{t\in \left[0,T\right]} $ independent, we obtain

$$ {\int}_0^T\varepsilon (t)\hskip0.1em dt\sim \mathcal{N}\left(0,T{\sigma}^2\right). $$

Thus $ {F}^{\star }(y) $ follows a Gaussian distribution centered on $ F(y) $ and with variance equals to $ T{\sigma}^2 $ . This makes it possible to compute confidence bounds on $ {F}^{\star }(y) $ . For a confidence level $ 1-u $ , $ u\in \left[0,1\right] $ , a confidence interval for $ {F}^{\star }(y) $ is obtained as

$$ {\mathtt{CI}}^{1-u}\left({F}^{\star }(y)\right)=F(y)\pm {\zeta}_{1-\frac{u}{2}}\sqrt{T}\sigma, $$

where $ {\zeta}_{1-\frac{u}{2}} $ is the quantile of order $ 1-\frac{u}{2} $ of the standard Gaussian distribution.

The assumption that $ f $ and $ {\sigma}^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.

4.1. Modeling

Here, the trajectories are supposed to be in a vertical plane and are defined by the altitude $ h $ , the Mach number $ \mathrm{M} $ and the engines rotational speed $ \mathrm{N}1 $ (expressed as a percentage of a maximal value). Hence, a trajectory $ y $ in this setting is a continuous $ {\unicode{x211D}}^3 $ -valued map defined on $ \left[0,T\right] $ , where $ T $ is a maximal climb duration fixed by the user. Hence, we have

$$ \forall t\in \left[0,T\right]\hskip2em y(t):= \left(h(t),\mathrm{M}(t),\mathrm{N}1(t)\right). $$

The quantity to minimize is the total fuel consumption $ \mathrm{TFC}:\mathcal{C}\left(\left[0,T\right],{\unicode{x211D}}^3\right)\to {\unicode{x211D}}_{+} $ which is defined via the fuel flow $ \mathrm{FF}:{\unicode{x211D}}^3\to {\unicode{x211D}}_{+} $ as followsFootnote ⁴:

$$ \mathrm{TFC}(y):= {\int}_0^T\mathrm{FF}\left(y(t)\right)\hskip0.1em dt. $$

Regarding the endpoints conditions, we require the trajectory to start at the altitude $ {h}_0 $ with Mach number $ {\mathrm{M}}_0 $ and to end at the altitude $ {h}_T $ with Mach number $ {\mathrm{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:

• • The rate of climb, that is the time-derivative of the altitude, has to be upper bounded by a given maximal value $ {\gamma}_{max} $ during the whole climb;

• • The Mach number should not exceed a certain value called the maximum operational Mach ( $ \mathrm{MMO} $ ).

The final time of the climb is given by $ {T}^{\star}\in \left[0,T\right] $ which is the first time where the aircraft reaches $ {h}_T $ with Mach number $ {\mathrm{M}}_T $ .

Finally, we mention that the fuel flow model $ \mathrm{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 Dewez et al. (Reference Dewez, Guedj and Vandewalle2020) 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.

4.2. 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 cannot be publicly released for commercial reasons. The data is here recorded by the quick access recorder (QAR).

Before considering the optimization 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 s to reduce the dataset 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\times {10}^{-2} $  kg/s and 1.73%.

Regarding the optimization, we are interested in climb phases from 3,000 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 optimization 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 $ \mathrm{N}1 $ values are let unconstrained. Finally the $ \mathrm{MMO} $ and $ {\gamma}_{max} $ are respectively set to 0.82 and 3,600 ft/min.

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 optimization problem. On the other hand, we use only the five 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 five 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 optimization problem in a space of dimension 20.

We are then in position to apply the optimization method developed in Section 2 using the PyRotor library. First of all a relevant value for $ {\nu}_{max}>0 $ has to be fixed. In order to propose a realistic optimized climb, we choose a $ {\nu}_{max} $ relatively small so that the optimized climb remains close to the reference ones. In particular, the quadratic objective function in (19) turns out to be convex for all $ \nu \in \left[0,{\nu}_{max}\right] $ permitting to use the quadratic programming solver from CVXOPT software imported in PyRotor. The preprocessing of the reference trajectories and the optimization 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 optimized trajectory obtained using PyRotor is given in Figure 3. We observe that the optimized 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 and Menéndez, Reference Codina and Menéndez2014 and references therein). In particular, the duration $ {T}^{\star } $ is equal to 1,033 s which is actually slightly shorter than the reference durations. We note also that the optimized Mach number shares a very similar pattern with the references. On the other hand, the optimized 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 optimized trajectory presents a realistic pattern inherited from the reference trajectories.

Figure 3. Optimized and reference altitudes, Mach numbers and engines rotational speeds—the optimized trajectory is represented by the blue curves.

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 optimized 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 modelings. In particular, the minimal percentage indicates the savings compared with the best trajectory (i.e., with the lowest consumption).

Table 1. Statistical description of the fuel savings of the optimized trajectory.

The savings are compared with the 48 recorded flights satisfying the present endpoints and the total consumption of the optimized trajectory is estimated sing the statistical model for the fuel flow. Q ₁, Q ₂, and Q ₃ refer to the first, second and third quartiles.

5.1. Modeling

To model this problem, we suppose without loss of generality that the trajectories are defined on the (time-) interval $ \left[0,1\right] $ and we let $ V:{\unicode{x211D}}^D\to {\unicode{x211D}}^D $ denote a vector field. Furthermore, the trajectories are assumed here to be continuously differentiable, that is they belong to $ {\mathcal{C}}^1\left(\left[0,1\right],{\unicode{x211D}}^D\right) $ . The work of $ V $ along a trajectory $ y\in {\mathcal{C}}^1\left(\left[0,1\right],{\unicode{x211D}}^D\right) $ is

$$ W\left(y,\dot{y}\right):= {\int}_0^1V{\left(y(t)\right)}^T\dot{y}(t)\hskip0.1em dt; $$

here $ \dot{y} $ 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\left(\dot{y}\right)={\int}_0^1\parallel \dot{y}(t){\parallel}_2^2\hskip0.1em 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 minimizing the cost function

$$ {F}_{\alpha}\left(y,\dot{y}\right):= \alpha J\left(\dot{y}\right)-W\left(y,\dot{y}\right)={\int}_0^1\alpha \parallel \dot{y}(t){\parallel}_2^2-V{\left(y(t)\right)}^T\dot{y}(t)\hskip0.1em dt, $$

where $ \alpha \geqslant 0 $ is arbitrarily chosen, is expected to lead to an optimized trajectory reflecting a trade-off between maximizing the work and minimizing the distance. Further we require the trajectory to stay in the hypercube $ {\left[0,1\right]}^D $ and to start and to end respectively at $ {y}_0\in {\left[0,1\right]}^D $ and $ {y}_1\in {\left[0,1\right]}^D $ .

Now, we remark that the above cost function involves the (time-)derivative $ \dot{y} $ . So one has to derive a formula permitting to compute the derivative of any trajectory $ {\left.y=\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\in {\mathcal{Y}}_{\mathcal{K}} $ from its associated vector $ c\in {\unicode{x211D}}^K $ , especially to compute $ \overset{\check{}}{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\in {\mathcal{Y}}_{\mathcal{K}} $ :

$$ \forall d=1,\dots, D\hskip2em {\dot{y}}^{(d)}=\sum \limits_{k=1}^{K_d}{c}_k^{(d)}{\dot{\varphi}}_k={\left({\left.\frac{d}{dt}\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\right)}^{(d)}. $$

We deduce then the following formula for $ \overset{\check{}}{F}(c) $ in the present setting:

$$ \overset{\check{}}{F}(c):= {F}_{\alpha}\left(\Phi {\left|{}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c,\frac{d}{dt}\Phi \right|}_{{\mathcal{Y}}_{\mathcal{K}}}^{-1}c\right). $$

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 $ \overset{\check{}}{F} $ explicit with respect to $ c $ in certain settings. For instance it is possible to derive an explicit quadratic formula for $ \overset{\check{}}{F}(c) $ when the integrand defining $ {F}_{\alpha } $ is quadratic with respect to $ y(t) $ and $ \dot{y}(t) $ ; this formula is implemented in PyRotor and the arguments to obtain it are similar to those proving Proposition 5.

5.2. Numerical results

Numerical results based on randomly generated data for the above physical application are presented in this section. We first consider trajectories with two components $ {y}^{(1)} $ and $ {y}^{(2)} $ lying in the square $ {\left[0,1\right]}^2 $ for the sake of simplicity. We set the starting and ending points as follows:

$$ {y}^{(1)}(0)=0.111\hskip1em ,\hskip1em {y}^{(2)}(0)=0.926\hskip1em ,\hskip1em {y}^{(1)}(1)=0.912\hskip1em ,\hskip1em {y}^{(2)}(1)=0.211 $$

with a tolerated error $ 1\times {10}^{-4} $ , and the vector field $ V:{\unicode{x211D}}^2\to {\unicode{x211D}}^2 $ is here defined by

$$ V\left({x}^{(1)},{x}^{(2)}\right)={\left(0,{x}^{(1)}\right)}^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 integrand of the cost function $ {F}_{\alpha } $ in the present setting is actually quadratic with respect to $ y(t) $ and $ \dot{y}(t) $ , so that an explicit quadratic formula for $ \overset{\check{}}{F}(c) $ implemented in PyRotor is available.

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 $ {\left[0,1\right]}^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}_{\alpha } $ .

In the present example, we set a $ {\nu}_{max} $ relatively large to explore a large domain around the reference trajectories. In this case, the objective function of the optimization problem (19) may be not convex even if it is still quadratic. So we make use of the generic optimization solver minimize(method=’trust-constr’) imported in PyRotor. Regarding the execution time, we have randomly and uniformly generated 100 values in the interval $ \left[0,10\right] $ for the parameter $ \alpha $ 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 Figure 4, we plot four optimized trajectories associated with different values of $ \alpha $ : 0, 0.35, 1, and 10. As expected the trajectory associated with the largest value of $ \alpha $ gives the most straight trajectory while the most curvy one is associated with $ \alpha =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 optimization approach may lead to optimized trajectories which slightly differ from the reference ones to reduce more the cost.

Figure 4. Optimized trajectories in the square $ {\left[0,1\right]}^2 $ for $ \alpha \in \left\{\mathrm{0,0.35,1,10}\right\} $ . Optimized 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 optimized trajectories and the bar shows the scale. Red arrows represent the pattern of the vector field $ V $ .

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

Table 2. Statistical description of the work gains in percentage for $ \alpha \in \left\{\mathrm{0,0.35,1,10}\right\} $

The values have been computed 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.