A Time-Dependent Subgraph-Capacity Model for Multiple Shortest Paths and Application to $\text {CO}_{\text {2}}$/Contrail-Safe Aircraft Trajectories

Abstract

This paper proposes a study motivated by the problem of minimizing the environmental impact of air transport considering the complete air network, thereby several aircraft. Both CO$_2$ and non-CO$_2$ effects are taken into account to calculate this impact. The proposed methodology takes a network point of view in which airspace capacities evolve as well as the traffic itself over time. Finding the shortest path with numerous constraints and various cost functions is a common problem in operations research. This study deals with the special case of multiple shortest paths with capacity constraints on a time-dependent subgraph. Multiple shortest paths are understood as one shortest path per vehicle considered. The static special case is modeled as a mixed integer linear program so that it can be efficiently solved by standard off-the-shelf optimization solvers. The time-dependent nature of the problem is then modeled via a sliding-window approach. Encouraging numerical results on the contrail-avoidance application show that the environmental impact can be significantly reduced while maintaining safety by satisfying the airspace capacity constraints.

Appendices

Appendix A: Time-Discretized Optimization Model

The mathematical optimization model obtained in the static case can be discretized to obtain the time-dependent optimization model. This transformation is based on the classical time-dependent shortest-path problem [45].

Some notations should be defined first:

• $T_i=\{ t_{0,i},\dots ,t_{f,i}\}$: the set of time slots for vehicle i;

• $T=\bigcup _{i=1}^{M} T_i$;

• $C_{k,t}$: the capacity of sector k, $k=1,2,\dots ,N$, for time slot t, $t\in T$;

• $w_{u,v,i}^{t}$: the cost for vehicle i, $i=1,2,\dots ,M$, to go through arc $(u,v)\in A$ at time slot t, $t\in T_i$;

• $\Delta _{u,v,i}^{t}$: the time necessary for vehicle i, $i=1,2,\dots ,M$, to go through arc $(u,v)\in A$ at time slot t, $t\in T_i$.

The decision variables are for each vehicle i, $i=1,2,\dots ,M$:

• $x_{u,v,i}\in \{ 0,1 \}$ is equal to 1 if vehicle i goes through arc (u, v);

• $z_{u,v,i}^{t}\in \{ 0,1 \}$ is equal to 1 if vehicle i enters arc (u, v) at time slot t;

• $y_{k,i}^{t}\in \{ 0,1 \}$ is equal to 1 if vehicle i flies through sector k at time slot t, $t\in T$.

The model is then:

$$\begin{aligned}\underset{X,Y,Z}{\min } & \sum _{i=1}^{M} \sum _{t\in T_i} \sum _{(u,v) \in A} w_{u,v,i}^{t} z_{u,v,i}^{t}&\end{aligned}$$

(A1a)

$$\begin{aligned}&\text {s.t.}&\sum _{(u,v) \in A} x_{u,v,i} - \sum _{(v,u) \in A} x_{v,u,i}=0,&~u \in V \backslash \{ s_i,e_i \},~ i=1,2,\dots ,M \end{aligned}$$

(A1b)

$$\begin{aligned} & \sum _{(s_i,v) \in A} x_{s_i,v,i} - \sum _{(v,s_i) \in A} x_{v,s_i,i}=1 ,&~ i=1,2,\dots ,M\end{aligned}$$

(A1c)

$$\begin{aligned} & \sum _{(e_i,v) \in A} x_{e_i,v,i} - \sum _{(v,e_i) \in A} x_{v,e_i,i}=-1 ,&~ i=1,2,\dots ,M\end{aligned}$$

(A1d)

$$\begin{aligned} & \sum _{(u,v) \in A} z_{u,v,i}^{t} - \sum _{(v,u) \in A} z_{v,u,i}^{t+\Delta _{v,u,i}^{t}}=0,&~u \in V \backslash \{ s_i,e_i \},~ i=1,2,\dots ,M, ~ t \in T_i\end{aligned}$$

(A1e)

$$\begin{aligned} & \sum _{(s_i,v) \in A} z_{s_i,v,i}^{t_{0,i}}=1 ,&~ i=1,2,\dots ,M\end{aligned}$$

(A1f)

$$\begin{aligned} & \sum _{t \in T_i} z_{u,v,i}^{t}=x_{u,v,i},&~ (u,v) \in A,~ k=1,\dots ,N, ~ i=1,\dots ,M \end{aligned}$$

(A1g)

$$\begin{aligned} & \sum _{i=1}^{M} y_{k,i}^{t} \le C_{k,t},&~ k =1,2,\dots ,N,~t\in T\end{aligned}$$

(A1h)

$$\begin{aligned} & y_{k,i}^{t}=1 \text { if and only if } \sum _{(u,v) \in A_k}z_{u,v,i}^{t}\ge 1 ,&~ i=1,2,\dots ,M , ~ k =1,2,\dots ,N, ~ t\in T_i \end{aligned}$$

(A1i)

$$\begin{aligned} & y_{k,i}^{t}=0,&~ i=1,2,\dots ,M , ~ k =1,2,\dots ,N, ~ t\in T \backslash T_i \nonumber \\ \end{aligned}$$

(A1j)

$$\begin{aligned} & X_i \in \{ 0,1 \}^{|A|},&~ i=1,2,\dots ,M\end{aligned}$$

(A1k)

$$\begin{aligned} & Z_{i,t} \in \{ 0,1 \}^{|A|},&~ i=1,2,\dots ,M,~ t \in T_i\end{aligned}$$

(A1l)

$$\begin{aligned} & Y_{i,t} \in \{ 0,1 \}^{N},&~ i=1,2,\dots ,M,~ t \in T_i. \end{aligned}$$

(A1m)

Constraints (A1b), (A1c), and (A1d) are the flow conservation constraints. Constraints (A1e) and (A1f) enforce consistency of space-time flow conservation. Constraints (A1g) make the link between the time-dependent and the static decision variables. Constraints (A1h) are the time-discretized capacity constraints. Constraints (A1i) and (A1j) define the auxiliary variables $y_{k,i}^{t}$. The former can easily be linearized, as for the static model using Proposition 1.

1.1 A.1 Comparison with the Proposed Heuristic

The resolution of this model via CPLEX [40] has been compared to the proposed heuristic approach. Table 3 shows the different results, with various instance sizes, and two capacity scenarios: one retricting and the other non-restricting, with a time-discretization step equal to the size of the sliding window. In this case, the time window is small, because the discretization step must be smaller than the flight time on the shortest arc. The impact of the value of the parameter K is also studied. These results are presented in Table 3.

Table 3 Comparison of the results obtained by solving directly the time-discretized model and using the sliding window approach for three instance sizes and two capacity scenarios

One observes differences between the two methods in terms of computation time and objective function value. On the one hand, the heuristic approach yields significantly lower computational times than that obtained by directly solving the time-discretized model. On the other hand, despite differences in the objective function value results, the values obtained with the heuristic approach remain close to those obtained with the exact method. The gain in computation time is what interests us here with regard to the application case. Thus, the loss in objective function is reasonable compared to the computational gain. All the more so as the chosen heuristic enables us to directly take into account the data update at each time interval, whereas solving the discretized model directly requires us to completely redo the computations each time. This low computation time is required when solving the problem in quasi-real-time framework but also a few hours before because computations can be repeated several times, to compare several scenarios for example, and it is essential that computation time is not too long.

Fig. 15

2D grid used for wind data extraction from Windy [46]

Appendix B: Data Processing for Numerical Experiments

This section shows how data for computational experiments have been extracted and computed. Appendix B.1 deals with wind data, while Appendix B.2 focuses on contrail data.

1.1 B.1 Wind Data

The costs defined by (7) involve the computation of the flight time over each arc $(u,v) \in A$. To compute these costs, the wind on the arcs, and the distance between u and v, the two ends of the arcs, have to be known.

Let $\lambda _u$ and $\lambda _v$ be the latitude of vertices u and v respectively, and let $\phi _u$ and $\phi _v$ be their longitude. The distance between u and v is given by:

$$\begin{aligned} d_{u,v}&=R~c_{rad}~\; \text {in km},\end{aligned}$$

(B1)

$$\begin{aligned}&=60~c_{degrees}~\; \text {in NM}, \end{aligned}$$

(B2)

where $R=6,371$ km is the Earth radius, $c_{rad}$ and $c_{degrees}$ represent the following c values, expressed in radians and in degrees respectively:

$$\begin{aligned} c = \arccos \Big (\sin (\lambda _u) \,\sin (\lambda _v) +\cos (\lambda _u)\,\cos (\lambda _v)\,\cos (\phi _v-\phi _u)\Big ). \end{aligned}$$

(B3)

Then, the flight time for aircraft i, noted $t_{u,v,i}$, between points u and v is given by:

$$\begin{aligned} t_{u,v,i}=\frac{d_{u,v}}{GS_{u,v,i}}, \end{aligned}$$

(B4)

where $GS_{u,v,i}$ is the ground speed of aircraft i on arc (u, v). It can be computed via:

$$\begin{aligned} GS_{u,v,i} = V_{a_{i}}+W_{u,v}, \end{aligned}$$

(B5)

where $V_{a_{i}}$ is the airspeed of aircraft i (considered constant), and $W_{u,v}$ is the wind encountered on the arc (u, v).

Wind data have been extracted from the website Windy [46] on a square grid of size $0.2^{\circ }$ above France, as shown in Fig. 15.

To compute the wind on each node of the graph, a so-called Shepard interpolation [47] was used. More precisely, for each node P located in a 2D-square $P_1P_2P_3P_4$ of the data grid (see Fig. 16), the wind W(P) at P is calculated from the wind at $P_k,~k=1,2,3,4$, using the distance from P to each of these points (noted respectively $d_1$, $d_2$, $d_3$ and $d_4$) as follows:

$$\begin{aligned} W(P)=\frac{\sum _{i=1}^{4} W(P_i) *d_i^{-p}}{\sum _{i=1}^{4} d_i^{-p}}, \end{aligned}$$

(B6)

where $p>1$ is a user-defined parameter (set to $p=2$ in this study).

Fig. 16

Notations for estimating the wind at P via Shepard interpolation

Finally, the wind along an arc (u, v) is simply defined as the average of that at u and at v:

$$\begin{aligned} W_{(u,v)}=\frac{W(u)+W(v)}{2}. \end{aligned}$$

(B7)

Figure 17 shows data used for the examples presented in the result section (Section 5).

1.2 B.2 Contrail Data

Contrails are formed in cold and humid areas. They persist and induce cirrus if the air is supersaturated in ice. The computation of persistent contrail areas is performed in two phases:

• 1. Areas favorable to contrail formation

• 2. Areas in which contrails will persist (ice supersaturated areas).

In [36], contrail areas are computed thanks to the Schmidt-Appleman criterion. This criterion gives a minimum threshold, $r_{min}$, of relative humidity of the air in liquid water, noted $RH_w$, above which contrails are formed: contrails are assumed to form when $RH_w \ge r_{min}$, where

$$\begin{aligned} r_{min}=\frac{G(T-T_{c})+e^{liq}_{sat}(T_{c})}{e^{liq}_{sat}(T)}, \end{aligned}$$

(B8)

$e^{liq}_{sat}(T)$ is the saturation vapor pressure over water, and $T_{c}$ is the estimated threshold temperature (in Celsius degrees) for contrail formation at liquid saturation. The latter is computed via:

$$\begin{aligned} T_{c} = -46.46+9.43\log (G-0.053)+0.72\log ^2(G-0.053), \end{aligned}$$

(B9)

where $G=\frac{EI_{H_2O}C_pP}{\epsilon Q(1-\eta )}$, $EI_{H_2O}=1.25$ is the water vapor emission index, $C_p=1004 ~J.kg^{-1}.K^{-1}$ is the heat capacity of the air, P is the ambient pressure (in Pascals), $\epsilon =0.6222$ is the ratio of the molecular masses of water and dry air, $Q=43*10^6 J.kg^{-1}$ is the specific heat of combustion, and $\eta = 0.3$ is the average propulsion efficiency of a commercial aircraft.

Fig. 17

Wind encountered in the example of the result

In [36], the ice super saturated areas are determined thanks to the following criterion: $RH_i>1$, where the relative humidity over the ice, noted $RH_i$ is computed as follows:

$$\begin{aligned} RH_i = RH_w * \frac{6.0612*\exp (\frac{18.102*T}{249.52+T})}{6.1162*\exp (\frac{22.577*T}{273.78+T})}, \end{aligned}$$

(B10)

and where T is the ambient temperature in Celsius degrees.

Fig. 18

Persistent-contrail areas in western France used for our instances

The relative humidity and temperature are also computed from data extracted from Windy [46] on a 2D-grid, and interpolated via quadratic interpolation. Figure 18 shows the data used for the examples presented in the result section (Section 5), where red areas are persistent-contrail-favorable areas, to be avoided.