2.2.1 Segment flight time prediction model

Given that the cruise phase constitutes most of the flight trajectory, focusing exclusively on optimising the cruise segment is reasonable. Assuming that the aircraft maintains a constant cruising altitude, let the coordinates of the start and end waypoints of segment ${a_i}$ are $\left( {{\lambda _i},{\phi _i}} \right)$ and $\left( {{\lambda _{i + 1}},{\phi _{i + 1}}} \right)$ , respectively. In this study, the distance and bearing between two waypoints are computed using the rhumb-line (loxodrome) formulation, which assumes a constant bearing between consecutive waypoints. The bearing angle ${\theta _i}$ and length ${l_i}$ of segment ${a_i}$ (see Fig. 1(c)) are defined as follows:

(2) \begin{align} \tan {\theta _i} = \frac{{{\phi _i} - {\phi _{i + 1}}}} {{\ln \left[ {\dfrac{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\lambda _i}}}{2}} \right)}}{{\tan \left( {\dfrac{\pi }{4} - \dfrac{{{\lambda _{i + 1}}}}{2}} \right)}}} \right]}} \end{align}

(3) \begin{align}{l_i} = \left\{ {\begin{array}{*{20}{c}}{\dfrac{{\left( {{R_E} + h} \right)\left( {{\lambda _{i + 1}} - {\lambda _i}} \right)}}{{\cos {\theta _i}}},\,} {}{{\lambda _{i + 1}} \ne {\lambda _i}}\\[3pt]{\left( {{R_E} + h} \right)\cos {\lambda _i}\left| {{\phi _{i + 1}} - {\phi _i}} \right|,\,} {}{{\lambda _{i + 1}} = {\lambda _i}}\end{array}} \right. \\[2 pt] \nonumber \end{align}

The wind vector along the segment ${a_i}$ can be decomposed into a parallel wind component ( ${\boldsymbol{\mathcal{W}}_{p,i}}\left( t \right)$ ) and a crosswind component ( ${\boldsymbol{\mathcal{W}}_{c,i}}\left( t \right)$ ). The crosswind is converted into an equivalent headwind to account for its impact on the effective ground speed. Accordingly, the ground speed ( ${V_{f,g,i}}\left( t \right)$ ) and the corresponding flight time ( ${t_{f,i}}$ ) of aircraft $f$ are calculated as follows:

(4) \begin{align}{V_{f,g,i}}( t ) = \sqrt {V_{f,T}^2 - W_{c,i}^2\left( t \right)} + {\boldsymbol{\mathcal{W}}_{p,i}}\left( t \right) \\[-28pt] \nonumber \end{align}

(5) \begin{align}\frac{{d{l_i}}}{{d{t_{f,i}}\left( t \right)}} = {V_{f,g,i}}\left( t \right) \\[5pt] \nonumber \end{align}

where, ${V_{f,T}}$ is the true airspeed. According to Equation (5), for each time instant $t \in \left[ {0,\boldsymbol{\mathcal{T}}} \right]$ , and for each ensemble member $\;k$ in the EPS, the corresponding flight time on segment ${a_i}$ can be represented as $t_{f,i}^k\left( t \right)$ . Consequently, the flight time ensemble is defined as: ${T_{f,i}}\left( t \right) = \left( {t_{f,i}^1\left( t \right),t_{f,i}^2\left( t \right), \cdots ,t_{f,i}^k\left( t \right), \cdots ,t_{f,i}^K\left( t \right)} \right)$ , where $K$ is the total number of ensemble members.

As previously noted, $\boldsymbol{\mathcal{W}}\left( t \right)$ is modelled as a piecewise function with discontinuities at time points $\left\{ {{t_v},2{t_v},3{t_v}, \cdots } \right\}$ . In this context, it becomes necessary to analyse the impact of such time-variant wind on consecutive flight segments. For example, suppose there are two consecutive segments ${a_i}$ and ${a_{i + 1}}$ , where aircraft $f$ is subject to distinct wind fields $\boldsymbol{\mathcal{W}}1$ and $\boldsymbol{\mathcal{W}}2$ , valid over the intervals $\left[ {0,{t_v}} \right]$ and $\left( {{t_v},2{t_v}} \right)$ , respectively. Let the expected arrival times of aircraft $f$ at the origin and destination of segments ${a_i}$ and ${a_{i + 1}}$ be denoted by $t_{f,i}^o$ , $t_{f,i}^d$ , $t_{f,i + 1}^o$ and $t_{f,i + 1}^d$ , respectively. Note that these segments are consecutive, so $t_{f,i}^d = t_{f,i + 1}^o$ . To determine the wind conditions experienced by aircraft $f$ on each segment, we use the expected arrival time at the origin of the segment as the basis. The following conditions are then established:

Condition 1: if $t_{f,i}^o \lt {t_v}$ and $t_{f,i}^d \lt {t_v}$ , both segments ${a_i}$ and ${a_{i + 1}}$ are influenced by the wind field $\boldsymbol{\mathcal{W}}1$ . Accordingly, ${T_{f,i}}\left( t \right) = T_{f,i}^1\left( t \right)$ and ${T_{f,i + 1}}\left( t \right) = T_{f,i + 1}^1\left( t \right)$ .

Condition 2: if $t_{f,i}^o \lt {t_v}$ and $t_{f,i}^d \ge {t_v}$ , then aircraft $f$ experiences a wind transition at the end of segment ${a_i}$ . In this case, ${T_{f,i}}\left( t \right) = T_{f,i}^1\left( t \right)$ and ${T_{f,i + 1}}\left( t \right) = T_{f,i + 1}^2\left( t \right)$ .

This approximation is reasonable because: (1) the wind field is defined as a piecewise random function with independent intervals, so crossing a discontinuity naturally triggers a new wind realisation; (2) keeping each segment under a single wind condition avoids the need for further sub-division within a segment, thereby maintaining computational tractability; and (3) in practical trajectory planning, the temporal resolution of wind forecasts (e.g., at the hourly scale) is coarser than the duration of a single flight segment, so assuming segment-wise constancy of wind is a widely adopted and reasonable engineering simplification. A small-sample test (see Appendix A) shows that forcing segment splits at TVMs changes average flight time by at most ±0.3% and airspace-level conflicts by at most ±0.5%, confirming that the baseline treatment (no split) yields negligible bias for the chosen discretisation.

This study assumes that each segment’s flight time is independent. This assumption is widely used in trajectory planning under wind forecast uncertainty [Reference Franco Espín, Rivas Rivas and Valenzuela Romero4, Reference Franco Espín, Rivas Rivas and Valenzuela Romero5, Reference Kamo, Rosenow, Fricke and Soler10, Reference Legrand, Puechmorel, Delahaye and Zhu11] and remains reasonable in the time-variant setting considered here. Since $\overrightarrow {\boldsymbol{\mathcal{W}}\left( t \right)} $ is modelled as a piecewise random function with discrete transitions, the wind experienced on different segments can be considered conditionally independent, especially when aircraft traverse them at distinct TVMs. Moreover, this assumption allows the model to remain tractable while preserving sufficient fidelity for planning purposes.

To further verify this assumption, a small-sample statistical check was conducted using simulated trajectories. The Pearson correlation coefficients between consecutive segment flight times were calculated across 100 randomly selected flights. The average correlation coefficient was below 0.1, indicating a weak correlation between adjacent segments. This result confirms that the independence assumption between segment flight times is reasonable within the scope of this study. Therefore, the expected arrival times can be computed as follows:

(6) \begin{align}{\mu _{f,i}}\left( t \right) = \frac{1}{K}\mathop \sum \limits_{k = 1}^K t_{f,i}^k\left( t \right)\end{align}

(7) \begin{align}t_{f,i}^o = {\sigma _f} + \Delta T + \mathop \sum \limits_{\bar a \in \bar A} \mu _{f,i}^{\bar a}\left( t \right)\end{align}

(8) \begin{align}t_{f,i}^d = t_{f,i}^o + {\mu _{f,i}}\left( t \right)\end{align}

where ${\sigma _f}$ denotes the departure time of aircraft $f$ , $\Delta T$ is the time from departure to the entry point of the target airspace, and $\bar A$ represents the set of segments traversed by aircraft $f$ before segment ${a_i}$ . In this study, $\Delta T$ is assumed to be a fixed constant for all flights, reflecting a uniform transition time before entering the planned airspace. Under this assumption, the $t_{f,i}^o$ and $t_{f,i}^d$ are primarily influenced by the departure time ${\sigma _f}$ and the cumulative expected flight times of the preceding segments. Therefore, in the following formulations, we use ${\sigma _f}$ as the representative variable to capture the temporal dynamics of each trajectory.

2.2.2 Time-variant shortest path problem

A lateral trajectory ( ${r_f}$ ) for aircraft $f$ can be represented as an ordered sequence of $n$ segments and their corresponding flight time ensembles, each of which depends on the departure time ${\sigma _f}$ , i.e.,

(9) \begin{align}{r_f} = \left\{ {\left( {{a_{f,1}},{T_{f,1}}\left( {{\sigma _f}} \right)} \right),\left( {{a_{f,2}},{T_{f,2}}\left( {{\sigma _f}} \right)} \right), \cdots ,\left( {{a_{f,i}},{T_{f,i}}\left( {{\sigma _f}} \right)} \right), \cdots ,\left( {{a_{f,n}},{T_{f,n}}\left( {{\sigma _f}} \right)} \right)} \right\},\\[2pt] \nonumber {a_{f,i}} \in A,\;{a_{f,i - 1}} \in A_{{a_{f,i}}}^ - ,\;\;{a_{f,i + 1}} \in A_{{a_{f,i}}}^ + ,\;{a_{f,i}} = \left( {{w_{f,i}},{w_{f,i + 1}}} \right),\;{w_{f,i}} \in W\end{align}

where, $A_{{a_{f,i}}}^ - $ denotes the set of segments entering waypoint ${w_{f,i}}$ , and $A_{{a_{f,i}}}^ + $ denotes the set of segments originating from waypoint ${w_{f,i + 1}}$ . Accordingly, the set of all segments associated with ${a_{f,i}}$ is defined as ${A_{{a_{f,i}}}} = A_{{a_{f,i}}}^ + \cup A_{{a_{f,i}}}^ - $ .

The lateral trajectory planning problem is formulated as a time-variant shortest path problem, aiming to minimise the expected total flight time under wind forecast uncertainty, subject to constraints on flight time deviation ${\delta _{f,i}}\left( {{\sigma _f}} \right)$ . The ${\delta _{f,i}}\left( {{\sigma _f}} \right)$ is defined as,

(10) \begin{align}{\delta _{f,i}}\left( {{\sigma _f}} \right) = \mathop {\max }\limits_{k \in K} \left\{ {t_{f,i}^k\left( {{\sigma _f}} \right)} \right\} - \mathop {\min }\limits_{k \in K} \left\{ {t_{f,i}^k\left( {{\sigma _f}} \right)} \right\}\end{align}

Based on graph theory, a time-variant directed graph $G\left( {W,A,{{\rm{\Phi }}_f}\left( {{\sigma _f}} \right),{{\rm{\Omega }}_f}\left( {{\sigma _f}} \right)} \right)$ can be constructed for a given aircraft $f$ at a fixed cruise altitude $h$ . In this graph,

• • $W$ represents the set of waypoints;

• • $A$ denotes the set of segments;

• • ${{\rm{\Phi }}_f}\left( {{\sigma _f}} \right) = \left\{ {{\mu _{f,i}}\left( {{\sigma _f}} \right)\;|\;{a_{f,i}} \in A} \right\}$ is the set of expected flight times on each segment under ${\sigma _f}$ ;

• • ${{\rm{\Omega }}_f}\left( {{\sigma _f}} \right) = \left\{ {{\delta _{f,i}}\left( {{\sigma _f}} \right)\;|\;{a_{f,i}} \in A} \right\}$ is the set of flight time deviations on each segment under ${\sigma _f}$ .

Then, the total expected flight time ${{\rm{\varphi }}_f}\left( {{\sigma _f}} \right)$ and the cumulative flight time deviation ${{{\unicode{x03C8} }}_f}\left( {{\sigma _f}} \right)$ of ${r_f}$ can be denoted as,

(11) \begin{align}{{\rm{\varphi }}_f}\left( {{r_f},{\sigma _f}} \right) = \mathop \sum \nolimits_{{a_{f,i}} \in {r_f}} {\mu _{f,i}}\left( {{\sigma _f}} \right),\,{{{\unicode{x03C8} }}_f}\left( {{r_f},{\sigma _f}} \right) = \mathop \sum \nolimits_{{a_{f,i}} \in {r_f}} {\delta _{f,i}}\left( {{\sigma _f}} \right)\end{align}

The optimal lateral trajectory $r_f^*$ that satisfies the constraint is,

(12) \begin{align}{\varphi _f}\left( {r_f^{\rm{*}},{\sigma _f}} \right) = \min \left[ {{\varphi _f}\left( {{r_f},{\sigma _f}} \right)} \right],\,\,{\rm{subject\, to}}\,{{{\unicode{x03C8} }}_f}\left( {{r_f},{\sigma _f}} \right) \le {{{\unicode{x03C8} }}_{f,0}}\end{align}

where, ${{{\unicode{x03C8} }}_{f,0}}$ is the predefined allowed deviation threshold.

2.3.1 Conflict detection model

Each aircraft must maintain a minimum separation of ${N_v} = 1,000\;{\rm{ft}}$ vertically and ${N_h} = 5NM$ laterally to ensure operational safety (see Fig. 2(a)) [Reference Nolan34]. A potential conflict is identified when the protection zones of two aircraft overlap, as illustrated in Fig. 2(b). It is important to note that such overlap does not necessarily indicate an actual conflict, but rather a situation that must be proactively avoided to maintain safety margins.

Figure 2. Aircraft separation criteria and potential conflict detection: (a) aircraft separation criteria; (b) potential conflict detection via zone overlap.

This study considers only temporal uncertainty and employs a four-dimensional (4D) grid-based conflict detection method to identify potential conflicts. To support this process, a 4D sampling point ensemble $\boldsymbol{\mathcal{S}_{\boldsymbol{f}}}\left({\boldsymbol{\sigma}}_{\boldsymbol{f}} \right)$ is generated based on the following principles: (1) the sampling points are uniformly distributed along the trajectory, with a constant spacing of $\Delta d$ between successive points; (2) the segment origins and destinations are included; (3) if no TVM exists on the segment, the sampling points are defined as shown in Fig. 3(a); (4) If a TVM is present, it must be explicitly included among the sampling points, as illustrated in Fig. 3(b). In Fig. 3, the red markers denote waypoints, the blue markers indicate TVM, and the yellow markers represent the generated sampling points.

(13) \begin{align}{\boldsymbol{\mathcal{S}}_{\boldsymbol{f}}}\left( {{\boldsymbol{\sigma}_{\boldsymbol{f}}}} \right) = \left\{ {{S_{f,1}}\left( {{\sigma _f}} \right),{S_{f,2}}\left( {{\sigma _f}} \right), \cdots ,{S_{f,s}}\left( {{\sigma _f}} \right), \cdots } \right\},\,{S_{f,s}}\left( {{\sigma _f}} \right) = \left( {{x_{f,s}},{y_{f,s}},{z_{f,s}},{T_{f,s}}\left( {{\sigma _f}} \right)} \right)\end{align}

Figure 3. The sampling principles according to the TVM: (a) segment without a TVM; (b) segment with a TVM.

where, ${x_{f,s}}$ , ${y_{f,s}}$ , and ${z_{f,s}}$ are the longitude, latitude, and altitude coordinates of sampling point $s$ , respectively; ${T_{f,s}}\left( {{\sigma _f}} \right)$ represents the predicted ensemble of times at which aircraft $f$ is expected to pass through the sampling point $s$ , given a departure time ${\sigma _f}$ .

To quantify the temporal dispersion at each sampling point, the maximum time deviation $D{T_{f,s}}$ is calculated as:

(14) \begin{align} &\qquad\qquad\qquad\quad D{T_{f,s}}\left( {{\sigma _f}} \right) = {\rm{max}}\left\{ {DT_{f,s}^{{\rm{min}}}\left({{\sigma _f}}\right)\!, DT_{f,s}^{{\rm{max}}}\left( {{\sigma _f}} \right)} \right\}\!, \\[2pt] \nonumber&DT_{f,s}^{{\rm{min}}}\left( {{\sigma _f}} \right) = {\mu _{f,s}}\left( {{\sigma _f}} \right) - \min \left\{ {{T_{f,s}}\left( {{\sigma _f}} \right)}\right\}\!,\,DT_{f,s}^{{\rm{max}}}\left( {{\sigma _f}} \right) = \max \left\{ {{T_{f,s}}\left( {{\sigma _f}} \right)} \right\} - {\mu _{f,s}}\left( {{\sigma _f}} \right) \end{align}

where, ${\mu _{f,s}}\left( {{\sigma _f}} \right)$ denotes the expected arrival time at point $s$ , calculated as the average of the ensemble ${T_{f,s}}\left( {{\sigma _f}} \right)$ . The predicted time window, denoted as $R{T_{f,s}}\left( {{\sigma _f}} \right)$ , is then defined as:

(15) \begin{align}R{T_{f,s}}\left( {{\sigma _f}} \right) = \left[ {{\mu _{f,s}}\left( {{\sigma _f}} \right) - D{T_{f,s}}\left({{\sigma _f}}\right)\!,{\mu _{f,s}}\left( {{\sigma _f}} \right) + D{T_{f,s}}\left( {{\sigma _f}} \right)} \right]\end{align}

This interval will be used to characterise the uncertainty in arrival time under time-variant wind conditions (see Fig. 4). Then, a conflict exists between two aircrafts ${f_1}$ and ${f_2}$ if either of the sampling points ${S_{{f_1},{s_1}}}$ and ${S_{{f_2},{s_2}}}$ violates one of the following conditions:

(16a) \begin{align}\sqrt {{{\left( {{x_{{f_1},{s_1}}} - {x_{{f_2},{s_2}}}} \right)}^2} + {{\left( {{y_{{f_1},{s_1}}} - {y_{{f_2},{s_2}}}} \right)}^2}} \gt {N_h} \\[-28pt] \nonumber \end{align}

(16b) \begin{align}\left| {{h_{{f_1},{s_1}}} - {h_{{f_2},{s_2}}}} \right| \gt {N_v} \\[-28pt] \nonumber \end{align}

(16c) \begin{align}R{T_{{f_1},{s_1}}} \cap R{T_{{f_2},{s_2}}} = \emptyset \end{align}

Figure 4. The predicted time window at sampling point $s$ .

2.3.2 Conflict resolution model

This study resolves potential conflicts by adjusting flight departure times. Traditional approaches typically aim to minimise the total number of conflicts or delays. However, optimising only one of these objectives may result in large delays for certain flights, compromising overall fairness. To address this limitation, a multi-objective optimisation model is developed that simultaneously considers the conflict numbers, flight delays and fairness.

Let the decision vector be defined as: $\boldsymbol{\sigma}_{\boldsymbol{f}} = \left( {{\sigma _{{f_1}}},{\sigma _{{f_2}}}, \cdots ,{\sigma _{{f_n}}}} \right)$ , where ${\sigma _{{f_i}}}$ represents the scheduled departure time of flight ${f_i}$ . The first objective is to minimise the total number of conflicts:

(17) \begin{align}{\rm{min\;}}{Z_{CN}} = \mathop \sum \limits_{i = 1}^{n - 1} \mathop \sum \limits_{j = i + 1}^n C\left( {{S_{{f_i}}}( {{\sigma _{{f_i}}}} ),{\rm{\;}}{S_{{f_j}}}( {{\sigma _{{f_j}}}} )} \right)\end{align}

where, the function $C( {{S_{{f_i}}}( {{\sigma _{{f_i}}}} ),\;{S_{{f_j}}}( {{\sigma _{{f_j}}}} )} ) = 1$ if the sampling point ensembles of flights ${f_i}$ and ${f_j}$ violate the separation criteria defined in Equation (16), and 0 otherwise.

Solving a three-objective optimisation problem – separately considering the conflict number, total delay and fairness – is computationally expensive and may lead to an overly complex Pareto front. To address this issue, we propose a unified metric called equilibrium delay (ED), which effectively balances total delay minimisation with fairness.

Let the flight delay ensemble be defined as:

(18) \begin{align} \boldsymbol{\mathcal{D}} = \left\{ {{D_{{f_1}}},{D_{{f_2}}}, \cdots ,{D_{{f_n}}}} \right\}\end{align}

where ${D_{{f_i}}} = {\sigma _{{f_i}}} - {\sigma _{{f_{i,0}}}}$ , ${\sigma _{{f_i}}}$ is the optimised departure time, ${\sigma _{{f_{i,0}}}}$ is the original scheduled departure time of flight ${f_i}$ . Define the mean and deviation of $\boldsymbol{\mathcal{D}}$ are ${\mu _D}$ and ${\delta _D}$ , respectively. Then, the second objective is to minimise the equilibrium delay, defined as:

(19) \begin{eqnarray} & \min {Z_{ED}} = \frac{{{\delta _D}}}{{{D_{{\rm{max}}}} - {\mu _D}}}& \nonumber\\ & {\mu _D} = \mathop \sum \limits_{i = 1}^n {D_{{f_i}}},{\delta _D} = {\rm{max}}\left\{ {{D_{{f_i}}}} \right\} - {\rm{min}}\left\{ {{D_{{f_i}}}} \right\}&\end{eqnarray}

where, ${D_{{\rm{max}}}}$ is the maximum allowable delay. The following constraint is imposed to ensure feasibility:

(20) \begin{align}{\sigma _{{f_i}}} - {\sigma _{{f_{i,0}}}} \le {D_{{\rm{max}}}}\,\,\forall {f_i} \in \boldsymbol{\mathcal{F}}\end{align}

where $\boldsymbol{\mathcal{F}}$ denotes the set of all flights.

Despite its simple formulation, ${Z_{ED}}$ provides an effective normalised metric that integrates both the ${\mu _D}$ and the ${\delta _D}$ , thereby capturing the balance between efficiency and fairness among flights. In this study, the ${D_{{\rm{max}}}}$ is set to 30 min. This value is commonly used in pre-tactical planning as a threshold beyond which departure adjustments may trigger further operational actions such as re-routing or tactical controller intervention. However, ${D_{{\rm{max}}}}$ is not an intrinsic constant; it can be flexibly adjusted according to the requirements of different airspace systems or operational policies.

The fairness metric ${Z_{ED}}$ is bounded within $\left[ {0,1} \right]$ . A smaller value indicates shorter average delays (higher efficiency) and smaller delay deviations (greater fairness). The case ${Z_D} = 0$ corresponds to ideal synchronisation where all flights depart without delay, whereas larger values reflect growing inequality in delay allocation. The behaviour of ${Z_{ED}}$ in edge cases is also intuitive. When the mean delay ${\mu _D}$ is close to zero but the delay variation ${\delta _D}$ is large, the denominator $\left( {{D_{{\rm{max}}}} - {\mu _D}} \right)$ remains near ${D_{{\rm{max}}}}$ while the numerator stays high, resulting in a relatively large ${Z_{ED}}$ . This correctly signifies a situation of low overall delay but poor fairness – an appropriate and desirable response for a fairness-oriented metric.

For illustration, consider three flights with delays of 2, 4 and 6 min, respectively. The average delay is ${\mu _D} = 4$ min, and the deviation ${\delta _D} = 1.63\;$ min. If, through optimisation, the delays become 3, 4 and 5 min, then ${\mu _D}$ remains 4 min but ${\delta _D}$ decreases to 0.82 min. Consequently, ${Z_{ED}}$ is reduced, indicating that the optimisation improves fairness while maintaining efficiency. This example highlights how ${Z_{ED}}$ captures both delay minimisation and equitable scheduling.

4.3.1 Sensitivity analysis of ${{\bf{\unicode{x03C8} }}_{f,0}}$

To assess the capability of ${{{\unicode{x03C8} }}_{f,0}}$ in characterising the uncertainty in trajectory flight time, a sensitivity analysis is conducted using a flight from Hong Kong International Airport (ICAO: VHHH) to Tel Aviv International Airport (ICAO: LLBG) as a case study. The entry and exit coordinates of the flight within CWA are $\left( {105.2^\circ E,28^\circ N} \right)$ and $\left( {78.8^\circ E,41.6^\circ N} \right)$ , respectively. The optimal lateral trajectories for the outbound (VHHH–LLBG) and return (LLBG–VHHH) flights are shown in Fig. 11, and the flight times are listed in Table 2.

As shown in Table 2, the flight time of downwind flights (LLBG–VHHH) is shorter than that of upwind flights (VHHH–LLBG), primarily due to the prevailing subpolar westerlies over the airspace. The wind field presented in Fig. 11, which is averaged over four hours, reveals a region of high uncertainty spanning $31.5^\circ - 34.5^\circ N$ and $95^\circ - 103^\circ E$ . This wind pattern leads to several observable trajectory phenomena.

When a higher level of reliability is required (i.e., a smaller value of ${{{\unicode{x03C8} }}_{f,0}}$ ), the trajectories tend to circumvent regions with high uncertainty. Notably, the trajectories (solid lines) of VHHH–LLBG are more inclined to pass directly through the regions with high uncertainty, whereas the trajectories (dashed lines) of LLBG–VHHH are more likely to divert around them. This directional difference is primarily attributed to the time-variant nature of the wind field: the uncertainty within the region intensifies over time, such that it remains relatively weak when the VHHH–LLBG flights enter the area, but becomes significantly stronger by the time the LLBG–VHHH flights pass through. In addition, the solid lines tend to deviate northwards to avoid strong, unfavourable headwinds – particularly in the region between $81^\circ E$ and $89^\circ E$ indicated by the green solid line – while dashed lines are more likely to take advantage of favourable tailwinds to shorten time, as illustrated by the yellow dashed line between $80.2^\circ E$ and $86^\circ E$ . These results demonstrate that the proposed model is sensitive to the parameter ${{{\unicode{x03C8} }}_{f,0}}$ , and effectively captures the influence of time-variant wind forecast uncertainty, indicating its potential for application in lateral trajectory planning.

4.3.2 Efficiency analysis of ${Z_{ED}}$

In the lower-level model, this study aims to minimise ${Z_{ED}}$ , thereby reducing total delay while maintaining fairness among flights. To evaluate the effectiveness of this objective, two commonly used strategies are selected for comparative analysis: minimising the mean delay time ( ${\mu _D}$ ) and minimising the delay range ( ${\delta _D}$ ). Table 3 presents the average values of ${\mu _D}$ and ${\delta _D}$ obtained from 20 independent runs for each strategy, based on the corresponding solutions with the minimum number of conflicts. To further validate the robustness of the proposed fairness measure ${Z_{ED}}$ , the Gini coefficient ( $\boldsymbol{\mathcal{G}}$ ) is also evaluated using the existing optimisation results. This indictor quantifies the overall dispersion and relative inequality of flight delays, providing complementary insight into fairness performance. As constraint (20) remains identical across all experiments, the number of conflicts in the minimum-conflict solutions is the same for all three strategies, ensuring comparability of results.

Table 3. Comparison of delay performance under different optimisation strategies; results are averaged over 20 independent runs

Table 4. Extremes of ${\mu _D}$ and ${\delta _D}$ for each strategy; results are averaged over 20 independent runs

Table 3 shows that when the strategy of minimising ${\mu _D}$ is adopted, the expected value $E\left( {{\mu _D}} \right)$ is the lowest among the three strategies, indicating high efficiency in reducing average delay. However, this strategy also produces the highest values of ${\delta _D}$ and $\boldsymbol{\mathcal{G}}$ , suggesting poor fairness among flights. In contrast, the strategy of minimising ${\delta _D}$ yields the lowest ${\delta _D}$ and $\boldsymbol{\mathcal{G}}$ , representing the fairest outcome, but at the cost of increased average delay. Conversely, the strategy of minimising ${Z_{ED}}$ achieves a balanced result, with both ${\mu _D}$ , ${\delta _D}$ , and $\boldsymbol{\mathcal{G}}$ values close to the optimal levels of the other two strategies. These results confirm that the proposed ${Z_{ED}}$ -based objective effectively reduces total delay while ensuring an equitable distribution of delay across flights.

Figure 11. Optimised lateral trajectories for flights VHHH–LLBG (solid lines) and LLBG–VHHH (dashed lines) under different ${{{\unicode{x03C8} }}_{f,0}}$ : ${{{\unicode{x03C8} }}_{f,0}} = 3$ (green), ${{{\unicode{x03C8} }}_{f,0}} = 5$ (red), ${{{\unicode{x03C8} }}_{f,0}} = 10$ (yellow).

Table 4 gives the extreme values of ${\mu _D}$ and ${\delta _D}$ across 20 independent runs under each strategy. When minimising ${\mu _D}$ , the value of ${\mu _D}$ remains relatively stable, but ${\delta _D}$ varies considerably, indicating difficulty in consistently ensuring fairness. Similarly, minimising ${\delta _D}$ results in greater variability in ${\mu _D}$ , reflecting a trade-off in efficiency. However, the ${Z_{ED}}$ -based strategy maintains both indicators within a narrow range, highlighting its robustness in achieving both delay reduction and fairness. These results further reinforce the practical value of the proposed objective function in multi-flight trajectory optimisation scenarios.

4.3.3 Efficiency analysis of SE strategy

To access the computational efficiency and convergence of the SE-based TVA^*, a comparative analysis was conducted against the distance-based TVA^*. Figure 12(a) presents the optimal trajectories obtained by both algorithms under ${{{\unicode{x03C8} }}_{f,0}} = 5$ . As illustrated, the green line generated by the distance-based TVA^* algorithm more closely follows the great circle path. However, despite its shorter distance, the corresponding flight time is longer, reaching 2.85 h. This is primarily attributed to the limited sensitivity of the distance-based TVA^* to uncertainty, which leads it to favour a shorter route while failing to account for the adverse effects of headwinds. In contrast, the SE-based TVA^* effectively incorporates time-variant wind forecast uncertainty and exhibits superior capability in mitigating the impact of unfavourable wind conditions. These results confirm the exceptional efficiency of the SE strategy.

Figure 12. Comparison between the SE-based and distance-based TVA^* algorithms, and the initial conflict distribution under ${{{\unicode{x03C8} }}_{f,0}} = 5$ : (a) optimal lateral trajectory of flight VHHH–LLBG: SE-based TVA^* (red) versus distance-based TVA^* (green); (b) initial optimal trajectories of all flights within the CWA, with potential conflicts indicated.

The flight time of the optimal trajectory calculated by GA is 2.72 h, which is identical to that produced by the SE-based TVA^*, thereby demonstrating the convergence of the proposed algorithm. Regarding computational performance, the GA requires 32 s to compute a single trajectory, whereas the SE-based TVA^* solves all 2,300 trajectories simultaneously in just 2.2 s, further validating its scalability and computational efficiency.

It is worth noting that the SE strategy is not strictly admissible or consistent under time-variant wind conditions, because the expected wind speed may slightly under- or over-estimate the true travel time of future segments when the wind field changes. Nevertheless, this does not affect the practical convergence behaviour of the algorithm. As demonstrated by the results above, the SE-based TVA^* attains exactly the same optimal trajectory as the GA benchmark, while achieving a substantial reduction in computation time. These results confirm that the lack of strict admissibility has negligible practical impact on the solution quality, and the SE-based TVA^* remains both reliable and highly efficient in the time-variant setting considered in this study.

Table 5. The average solution time of MOCCEA based on EG, RG, DG and DEG strategies (min); results are averaged over 20 independent runs

4.3.4 Efficiency analysis of DEG strategy

The effectiveness of the DEG strategy is compared with three representative existing grouping methods: equilibrium grouping (EG), random grouping (RG), and dynamic grouping (DG) [Reference Guan, Zhang, Lv, Chen and Michal16] across four different traffic scenarios. The number of flights in each scenario varies due to differences in the time span: 465 flights (08:00–09:00), 949 flights (08:00–10:00), 1,708 flights (08:00–11:00) and 2,300 flights (08:00–12:00), respectively. Figure 12(b) illustrates the distribution of the initial optimal trajectories and conflict locations for the 2,300-flight scenario. The red lines denote the trajectories, while the yellow markers indicate the detected conflicts. A total of 2,182 conflicts have been identified, predominantly concentrated in the eastern and northern regions of the airspace, where traffic density is highest.

Table 5 presents the average solution times for each strategy, based on 20 independent runs. The results indicate that DEG consistently outperforms the other methods, with the most pronounced advantage observed in large-scale problems. This performance gain can be attributed to the high degree of coupling among large flights, which complicates the balanced allocation of computational load across CPU threads, thereby limiting the full parallel efficiency of MOCCEA.

Figure 13. The efficiency and convergence analysis: (a) the size of the largest group and the number of groups during the optimisation process of the DG strategy where $n = 2,300$ ; (b) the IGD trends of the DEG-based MOCCEA across the four scenarios.

To further illustrate this issue, Fig. 13(a) presents the evolution of the maximum group size and the number of groups during the optimisation process of the DG strategy for the case where $n = 2,300$ . The DG strategy performs grouping based solely on pairwise flight conflict relationships, without considering computational load balancing. As the algorithm proceeds, the size of the largest group gradually decreases, while the number of groups tends to increase. This trend reflects the reduction in conflict density, as more flights become conflict-free and interdependencies are weakened. However, groups containing only a single flight persist throughout the process. Although they no longer contribute significantly to optimisation, they consume computational resources, reducing overall efficiency. The DEG strategy effectively addresses this limitation by incorporating conflict structure and load balancing into the grouping process. For the largest scenario with 2,300 flights, the DEG strategy reduces the computation time to 42.28 min, achieving speed-ups of 1.42× over EG (60.08 min), 1.23× over RG (52.03 min), and 1.16× over DG (49.22 min).

As discussed in Section 4.3.3, the SE-based TVA^* requires only 2.2 s to compute all 2,300 trajectories, accounting for less than 1% of the total solution time. The DEG-based MOCCEA dominates the overall computational cost owing to the large-scale conflict detection and resolution process. In serial execution, solving the entire bi-level optimisation problem takes approximately 266.4 min, whereas the parallel implementation using eight threads on an eight-core CPU reduces the total computation time to 42.28 min, achieving an average speed-up of 6.3×. These results clearly demonstrate that the proposed parallel framework substantially enhances computational efficiency and scalability, making it well-suited for large-scale strategic trajectory planning applications.

Figure 13(b) illustrates the IGD convergence trend of the DEG-based MOCCEA across the four scenarios, showing that the algorithm typically converges within approximately 120 generations. These findings confirm that the DEG-based MOCCEA is convergent and significantly improves computational efficiency, particularly in large-scale, high-density traffic scenarios.

Figure 14 presents the Pareto fronts for the four scenarios under comparison. Across all scenarios, the DEG strategy consistently produces superior Pareto fronts compared to the other strategies, demonstrating its overall advantage. Specifically, the DG strategy increases the number of groups during the conflict resolution process, which results in additional computational overhead. The EG strategy, by contrast, overlooks the coupling relationships among flights, leading to inefficient grouping and reduced performance. In comparison, the DEG strategy effectively balances computational loads across threads while preserving critical conflict structures, enhancing both the final solutions’ efficiency and quality.

To better emphasise the methodological contribution of the proposed DEG strategy, Table 6 presents a comparative analysis with the EG, RG and DG methods. Unlike EG and RG, which are static in nature, or DG, which groups flights solely according to conflict relationships, the proposed DEG method integrates both conflict awareness and thread load balancing. This dual consideration enables DEG to maintain near-uniform computational workloads across threads while preserving critical inter-flight conflict structures, thereby achieving superior scalability and parallel efficiency in large-scale optimisation.

Table 6. Comparison between DEG and existing grouping strategies

4.3.5 Evaluation of the bi-level STP framework

To demonstrate the advantages of the proposed bi-level planning model, this section compares the optimal Pareto fronts generated by the bi-level model and the two-stage model under four traffic scenarios, as illustrated in Fig. 15. The results show that, across all scenarios, the Pareto fronts obtained from the bi-level model consistently dominate those of the two-stage approach. This performance gain is primarily attributed to the inherent limitations of the two-stage model: under time-variant wind conditions, adjusting departure times during the second stage may render the lateral trajectories optimised for the first stage suboptimal. Moreover, changes in departure times reduce the predictability of aircraft arrival times at trajectory sampling points, making it difficult to accurately detect and resolve potential conflicts. In contrast, the bi-level model accounts for the coupling between trajectory planning and conflict resolution, thereby effectively addressing these issues.

Table 7. Average solution time (min) of the two-stage model and bi-level models across four traffic scenarios (8 CPU threads); results are averaged over 20 independent runs

Figure 14. The Pareto front of the EG, RG, DG and DEG strategies for four scenarios: (a) $n = 465$ ; (b) $n = 949$ ; (c) $n = 1,708$ ; (d) $n = 2,300$ .

Figure 15. The Pareto fronts of the two-stage model and bi-level model for four scenarios: (a) $n = 465$ ; (b) $n = 949$ ; (c) $n = 1,708$ ; (d) $n = 2,300$ .

Table 7 shows the computation times for both the two-stage and bi-level models across the four scenarios. As the number of flights increases, the bi-level model exhibits increasingly better performance compared with the two-stage model, highlighting its scalability and robustness in handling large-scale optimisation. These results confirm that the proposed approach is both practical and effective for strategic trajectory planning, with strong potential to enhance the overall efficiency of air traffic management.