2.1 Nonlinear continuous model

To simplify the model studied herein, Two-Degree-of-Freedom (2DOF) point-mass equations of motion with a non-rotating spherical Earth are considered in this paper. The ascent phase dynamics can thus be written as^{(Reference Miele25)}

(1) $ \begin{eqnarray}m\dot V = T\cos\alpha-D-mg\sin\gamma \end{eqnarray}$

(2) $ \begin{eqnarray}mV\dot \gamma =T\sin\alpha+L-mg\cos\gamma+mV^2/r\cos\gamma\end{eqnarray}$

(3) $ \begin{eqnarray}\dot m = -\frac {{T}}{{I_{sp}.g_0}} \end{eqnarray}$

(4) $ \begin{eqnarray} \dot r = V\sin \gamma \end{eqnarray}$

(5) $ \begin{eqnarray}\dot \beta = \frac {{ V\cos \gamma}}{{r}}\end{eqnarray}$

where T, D and L are the thrust, drag and lift forces, respectively, $$ D = qSC_D $$ and $$ L = qSC_L $$ , where $$ q = 1/2\rho V^2 $$ is the dynamic pressure. The density $$ \rho $$ is modelled using the International Standard Atmosphere model. V is the velocity, $$ \gamma $$ is the fight path angle, $$ \beta $$ is the range angle, m is the sum of the masses of the launch vehicle and the payload, r is the distance between the launch vehicle and the Earth centre, g is the local gravitational acceleration and $$ I_{sp} $$ is the specific impulse.

The Earth is taken to be a non-rotating, spherical body whose surface is described by the mean sea-level radius and whose gravitational acceleration varies with radius r as follows:

(6) $ \begin{equation}g =g_0 \left(\frac{{R_0 }}{{r}}\right)^2\end{equation}$

where $$ R_0=6378.137\mathrm{km} $$ is the mean radius of the Earth.

The thrust model is given by

(7) $ \begin{equation}T=T_{vac}-A_e \cdot P,\end{equation}$

where $$ T_{vac} $$ is the engine thrust in vacuum, P is the atmospheric pressure at the current altitude of the SLV and $$ A_e $$ is the exit area of the engine nozzle.

The SLV is a two-stage solid rocket engine. The take-off mass of the SLV consists of the inert vehicle mass, the propellant mass, the payload mass, the payload margin mass and the payload fairing mass (Table 1).

Table 1

Mass and engine characteristics

The aerodynamic model of the vehicle predicts the total coefficients of lift $$ C_L $$ and drag $$ C_D $$ as functions of the angle-of-attack $$ \alpha $$ and Mach number. Using the input aerodynamic data relationship between $$ \alpha $$ and $$ C_L $$ ( $$ C_L=f(\alpha,Mach) $$ ) in combination with Equation (2), the trimmed angle-of-attack corresponding to the required lift coefficient can be obtained implicitly at any time (or Mach number).

2.2 Flight model constraints

The SLV performance is subject to physical constraints that affect its manoeuvrability. During the denser atmospheric flight of the vehicle, it is essential to ensure that the structural loads experienced by the vehicle remain well within the permissible limits. Here, the product of dynamic pressure and aerodynamic angle-of-attack has been constrained to be well within the maximum allowable value as expressed below:

(8) $ \begin{equation}\begin{array}{l} \left| {q\alpha} \right| \le {(q\alpha)_{max}} \\ \end{array}\end{equation}$

Also, there is an equation for the angle-of-attack $$ \alpha $$ as a nonlinear constraint. For the programmed flight, this is well described by using the following equation:

(9) $ \begin{equation}\begin{array}{l} \alpha(t) = \left\{ {\begin{array}{*{20}{c}}{0}&\quad {if}&\quad {t\le t_0}\\ \\[-7pt] 4\alpha_{\max}e^{-a(t-t_0)}(e^{-a(t-t_0)}-1)&\quad {if}&\quad {t_0<t\le t_1}\\ \\[-7pt]-4\alpha_{\max}e^{-a(t-t_1)}(e^{-a(t-t_1)}-1)&&\quad otherwise\\\end{array}} \right.\\\end{array}\end{equation}$

where $$ t_0 $$ is the vertical flight time, $$ t_1 $$ is the burn-out time of the rocket engine of stage 1, $$ \alpha_{\max} $$ is the maximum angle-of-attack and a is a positive constant. Moreover, $$ T_a=1/a $$ is the time constant, which should be chosen depending on the type of SLV to achieve the maximum value of the angle-of-attack at a desired time of flight or a zero angle-of-attack before the specified time (before transonic flight).

Note that the pitch angle may be computed as

(10) $ \begin{equation}\theta(t)=\alpha(t)+\gamma(t).\end{equation}$

2.3 The cost functional

For trajectory optimisation in multi-stage launch problems, a typical cost objective is to maximise the payload mass. In this work, a weighted sum functional with components of the final states of the velocity V and the flight path angle $$ \gamma $$ as well as an integral term, namely minimised flight time (or maximised final mass), is considered:

(11) $ \begin{equation}{C}\left( x \right) = \sum\limits_1^{v} {{W_i}{C^i(x)}} ={W_V}{C^{V}} + {W_\gamma}{C^{\gamma}}+ {W_t}{C^{t}}\\ = {W_V}(V-V_f) + {W_{\gamma}}(\gamma-\gamma_f)+{W_t}\int\limits_{{t_0}}^t {dt} \end{equation}$

where v is the number of objective functions, p is the number of inequality constraints and $$ x\in R^n $$ is a vector of design variables, where n is the number of independent variables $$ x_v $$ . $$ C(x)\in R^v $$ is a vector of objective functions $$ C^i(x):R^n\to R^1 $$ , $$ i=1,..,v $$ . $$ W_t, W_V $$ and $$ W_{\gamma} $$ represent the weights of the flight time, velocity and flight path angle penalty criterion, respectively. These weights are typically set by the decision-maker such that they sum to unity.

Moreover, in this optimisation study, $$ t_0 $$ and $$ \theta_S $$ are chosen as design variables, i.e. $$ x=[t_0;\theta_S] $$ .

Note that, if different objective functions have different units, the norm of any degree becomes insufficient to achieve mathematical closure. Consequently, the objective functions should be transformed such that they become dimensionless. A common function transformation used for this is presented below^{(Reference Marler and Arora22)}:

(12) $ \begin{equation}{\bar C^i} = \frac{{{C^i}}}{{C^i_{\max }}}\end{equation}$

which results in a non-dimensional objective function with an upper limit of 1 (note that $$ C^i_{\max } \ne 0 $$ is assumed). The detailed computation of $$ C^i_{\max } $$ is discussed in the next section.

In the optimal trajectory planning, minimising the cost functional in Equation (11) is of interest.

3.1 Optimal discrete model definition

In practice, the model over a special region is discretised, which includes the spatial position, and the trajectory optimisation problem is formulated as a constrained MCNF problem over a layered network flow structure. In this work, a digitised grid around the Earth is used to represent the path planning environment. This grid network with an $$ m\times n $$ grid is used as an input to the algorithm. Each grid point is of the polar form $$ (r,\beta) $$ , representing the position coordinates with respect to the Earth-centred reference frame. A network flow structure is constructed based on the grids. That is, each node $$ (i,\,j) $$ of the network represents the position coordinates of the form $$ (r_{(i,\,j)},\beta_{(i,\,j)}) $$ . The nodes can be classified as a layered network with n layers, with each layer containing m nodes. The first layer is at ground level (launch altitude), while the last layer is at orbital height. Naturally, there are $$ m\times n $$ nodes. Since the aim is to find the optimised path, the SLV should always fly toward the goal node. Thus, each node in the $$ j^{th} $$ layer is connected to the K nodes in the $$ (\,j-1)^{th} $$ layer via K arcs, where K represents the number of nodes in a predefined field of view, so arc $$ [(k,\,j-1),(i,\,j)] , k=1,2,...,K, K\le m $$ connects the $$ k^{th} $$ node in the $$ (\,j-1)^{th} $$ layer as the head of the arc to the $$ i^{th} $$ node in the $$ j^{th} $$ layer as the tail of the arc, as shown in Figure 1. For instance, in the figure, we have $$ K=7 $$ .

Figure 1.

The layered network constructed around the Earth.

Consider a directed network $$ G = (N, A) $$ , where N is a set of $$ m\times n $$ nodes and A represents a set of $$ m\times n\times K $$ arcs. Note that all of the arcs connect nodes in two consecutive layers, as mentioned above. The MCNF problem is to find a feasible flow such that the sum of the transmission costs is minimum, which can be formulated as follows^{(Reference Ahuja, Magnanti and Orlin26)}:

(13) $ \begin{equation}\begin{array}{l} Minimise\,J = \sum\limits_{\{ (i,\,j) \in N\} }\,\, {\sum\limits_{\{ k \in K\} } {{C_{k,\,j - 1,i,\,j}}{X_{k,\,j - 1,i,\,j}}} } \\ \\[-7pt] S.T.\,\left\{ {\begin{array}{*{20}{c}}{\sum\limits_{\{ [(k,\,j - 1),(i,\,j)] \in A\} } {{X_{k,\,j - 1,i,\,j}}} - \sum\limits_{\{ [(k,\,j - 1),(i,\,j)] \in A\} } {{X_{i,\,j,k,\,j - 1}}} ={b_{i,\,j}}} \\ \\[-7pt] {0 \le {X_{k,\,j - 1,i,\,j}} \le {U_{k,\,j - 1,i,\,j}},\,\,\,\,\{ [(k,\,j - 1),(i,\,j)] \in A\} } \\\end{array}} \right. \\ \end{array}\end{equation}$

where $$ X_{k,\,j-1,i,\,j} $$ denotes the flow on arc $$ [(k,\,j-1),(i,\,j)] $$ , $$ C_{k,\,j-1,i,\,j} $$ is the cost to transport one unit of flow on that arc and $$ b_{i,\,j} $$ denotes the flow produced or consumed at node ( $$ i,\,j $$ ). If $$ b_{i,\,j}> $$ 0, node ( $$ i,\,j $$ ) is called a source (supply) node, while a node ( $$ i,\,j $$ ) with $$ b_{i,\,j}<0 $$ is called a sink (demand) node. Otherwise (i.e., $$ b_{i,\,j} = 0 $$ ), node ( $$ i,\,j $$ ) is called a transshipment node. $$ U_{(k,\,j-1,i,\,j)} $$ denotes the upper bound on the flow on each arc $$ [(k,\,j-1),(i,\,j)]\in A $$ .

The first constraint in Equation (13) is called the flow conservation or mass balance constraint. The mass balance constraint states that the outflow minus inflow must equal the supply/demand of the node. If the node is a supply node, its outflow exceeds its inflow; if the node is a demand node, its inflow exceeds its outflow; and if the node is a transshipment node, its outflow equals its inflow. The flow must also satisfy the lower bound and capacity constraints (the second constraint of Equation 13), which we refer to as flow bound constraints. The flow bounds typically model physical capacities or restrictions imposed on the flows’ operating ranges. A set of values for the flow variables X satisfying all the constraints is called a feasible point or a feasible solution. The set of all such points constitutes the feasible region or the feasible space.

3.1.2 The optimality conditions

To solve the optimisation problem in Equation (13), after determining the feasible space, the cost function must be minimised. A variant of the MCNF algorithm is proposed, which is suitable for layered networks. This algorithm uses the optimality conditions of a modified shortest-path problem to find the optimum path. To express the optimality conditions, a cost label is associated with each node of the network, which represents the amount of cost from the start node to that node. The cost label of node ( $$ i,\,j $$ ) is referred to as p( $$ i,\,j $$ ).

Optimum-path algorithms typically rely on the property that an optimum path between two nodes contains other optimum paths within it. This optimal substructure property is a hallmark of the applicability of the proposed algorithm. The following lemma states the optimal substructure property of optimal paths more precisely:

Lemma 1. Subpaths of optimum paths are optimum paths; Let $$ P = \left\langle {{v_1},{v_2},...,{v_n}} \right\rangle $$ be an optimum path from node $$ v_1 $$ to node $$ v_1 $$ and, for any l and t such that $$ 1 \le l \le t \le n $$ , let $$ {P_{lt}} = \left\langle {{v_l},{v_{l + 1}},...,{v_t}} \right\rangle $$ be the subpath of P from node $$ v_l $$ to node $$ v_t $$ . Then, $$ {P_{lt}} $$ is an optimum path from $$ v_l $$ to $$ v_t $$ .

Proof. Trivial.

Theorem 1. Using the following theorem, the cost optimality can be shown simply as well:

The necessary and sufficient condition for the optimality of the cost labels is that, for any node ( $$ i,\,j $$ ) connected to node $$ (k,\,j-1) $$ via arc $$ [(k,\,j-1),(i,\,j)] $$ ,

(14) $ \begin{equation}p(i,\,j) \le p(k,\,j - 1) + {C_{k,\,j-1,i,\,j}}.\end{equation}$

Moreover, arc $$ [(k,\,j-1),(i,\,j)] $$ is on the optimal path if and only if

(15) $ \begin{equation}p(i,\,j) = p(k,\,j - 1) + {C_{k,\,j-1,i,\,j}}.\end{equation}$

Proof. The interested reader can find the proof of the theorem in reference^{(Reference Zardashti23)}.

3.2 The proposed algorithm

In essence, the algorithm proceeds from the start (or initial) node. Then, at each step, nodes of one layer are tagged with their optimal cost label until the goal node is reached. As will be seen, the sequential structure of this scheme will allow us to introduce the algorithm.

3.2.1 Computing the optimal cost labels

To start the algorithm, the cost label of the start node is set to zero, then the cost label of any other nodes is set to infinity. Then a forward scheme is applied to calculate the cost labels of nodes in each layer, considering the optimality conditions and the cost labels of the previous layer. Supposing that the optimal cost labels of the $$ (\,j-1)^{th} $$ layer have been computed, for each node k in $$ (\,j-1)^{th} $$ layer, then the value of any constrained parameter described in Equation (11) is calculated on the arc $$ [(k,\,j-1),(i,\,j)] $$ . If the resulting values do not satisfy their bounds, the flow on the related arc is set to zero ( $$ {X_{k,\,j-1,i,\,j}} = 0 $$ ). If the flows on arcs connected to node ( $$ i,\,j $$ ) are zero, p( $$ i,\,j $$ ) remains at infinity, otherwise:

(16) $ \begin{equation}p(i,\,j) = \mathop {\min }\limits_{k:{X_{k,\,j-1,i,\,j}}\ne\ 0} [p(k,\,j - 1) + {C_{k,\,j-1,i,\,j}}].\end{equation}$

In this case, $$ {X_{k_{opt},j-1,i,\,j}} $$ is set to 1, where $$ k_{opt}=\mathop {\arg \min }\limits_{k:{X_{k,\,j-1,i,\,j}}\ne\ 0} [p(k,\,j - 1) +C_{k,\,j-1,i,\,j}] $$ . Moreover, node $$ (k_{opt},j-1) $$ represents the parent (predecessor) of node ( $$ i,\,j $$ ), that is, $$ (k_{opt},j - 1) = PARENT(i,\,j) $$ .

The expression $$ {X_{k_{opt},j-1,i,\,j}} = 1 $$ indicates that the solution to the problem will send 1 unit of flow from the starting node to any node along the optimum path.

According to Equations (14) and (15), p( $$ i,\,j $$ ) is the optimal cost label of node ( $$ i,\,j $$ ). Because this smallest label node will have its label permanently set and never revised again, we still refer to this procedure as a label-setting algorithm.

Following this scheme, the optimal cost label of all nodes including the goal node can be computed. Note that, since the computed cost labels satisfy the optimality conditions in Equations (14) and (15), they represent the optimum cost $$ C_{k,\,j-1,i,\,j} $$ from the start node to those nodes. Specially, the label of the goal node represents the total cost of the optimum path from the start node to the goal node.

3.3 Discrete model formulation

To implement the layered network using network search, a discrete space of the SLV model must be provided. The length (Euclidean distance), the flight path angle and the heading angle of the arc $$ [(k,\,j-1),(i,\,j)] $$ may be computed as follows:

(17) $ \begin{align} \Delta {\beta} & =\beta_{(i,\,j)}-\beta_{(k,\,j - 1)}\nonumber\\[5pt] \Delta {s_{(k,\,j - 1,i,\,j)}} & = \sqrt {r_{(i,\,j)}^2 + r_{(k,\,j - 1)}^2 -2{r_{(i,\,j)} r_{(k,\,j - 1)}} cos(\Delta {\beta})} \\[5pt] \gamma _{(k,\,j - 1,i,\,j)} & = pi/2- {\sin ^{ - 1}\Big[\frac{{r_{(i,\,j)}}} {\Delta {s_{(k,\,j - 1,i,\,j)}}}}\sin(\Delta {\beta})\Big]\nonumber\end{align}$

By replacing time versus length, one gets

(18) $ \begin{equation}{V_{(k,\,j - 1)}} = \frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{\Delta t}} \Rightarrow \Delta t = \frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j -1)}}}}\end{equation}$

An expression for $$ q\alpha $$ may be obtained as

(19) $ \begin{equation}q\alpha_{(k,\,j - 1)} =1/2\rho_{(k,\,j-1)} {V_{(k,\,j - 1)}^2}.{\alpha _{(k,\,j - 1)}}\end{equation}$

where $$ \alpha _{(k,\,j - 1)} $$ is obtained using the equation

(20) $ \begin{equation}\begin{array}{l} \alpha _{(k,\,j - 1)} = \left\{ {\begin{array}{*{20}{c}}{0}&{if}&{{t_{(k,\,j - 1)}}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_0}}\\\\[-7pt] 4\alpha_{max} e^{-a[t_{(k,\,j - 1)}-t_0]}(e^{-a[t_{(k,\,j - 1)}-t_0]}-1)&{if}&\quad{{t_0} < {t_{(k,\,j - 1)}}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_1}}\\\\[-7pt] -4\alpha_{max} e^{-a[t_{(k,\,j - 1)}-t_1]}(e^{-a[t_{(k,\,j - 1)}-t_1]}-1)&&otherwise\\\end{array}} \right.\\\end{array}\end{equation}$

and $$ \theta _{(k,\,j - 1)} $$ is computed as

(21) $ \begin{equation}\theta_{(k,\,j - 1)} =\gamma _{(k,\,j - 1)}+\alpha _{(k,\,j - 1,i,\,j)}.\end{equation}$

An expression for the rate of the flight path is

(22) $ \begin{equation}\begin{array}{l} {{\dot \gamma }_{(k,\,j - 1)}} \cong \dfrac{{\Delta {\gamma _{(k,\,j - 1,i,\,j)}}}}{{\Delta {s_{(k,\,j -1,i,\,j)}}}}{V_{(k,\,j - 1)}}, \end{array}\end{equation}$

where $$ \Delta {\gamma _{(k,\,j - 1,i,\,j)}} = {\gamma _{(k,\,j - 1,i,\,j)}} - {\gamma _{(kk,\,j - 2,k,\,j - 1)}} $$ .

The lift and drag forces are obtained (by using Equation 2) implicitly as

(23) $ \begin{equation}\begin{array}{l}{L_{(k,\,j - 1)}} = {{m_{(k,\,j - 1)}}.{V_{(k,\,j - 1)}}.{{\dot \gamma }_{(k,\,j - 1)}} + {m_{(k,\,j - 1)}}.g.\cos {\gamma _{\left( {k,\,j - 1}\right) \to \left( {i,\,j} \right)}}} - {T_{(k,\,j - 1)}}\sin {\alpha _{(k,\,j - 1)trim}} \\ \\[-7pt]{CD _{(k,\,j - 1)trim}} = f(C{L_{(k,\,j - 1)}}). \\ \end{array}\end{equation}$

The velocity and mass states are updated as follows:

(24) $ \begin{equation}\begin{array}{l}{T_{(k,\,j - 1)}} = \left\{ {\begin{array}{*{20}{c}}{{T_1}}&\quad{if}&{{t_{(k,\,j - 1) }}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_1}}\\ \\[-7pt]{{T_2}}&\quad {if}&{{t_1} < {t_{(k,\,j - 1) }}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_2}}\\ \\[-7pt]0&{}&{otherwise}\end{array}} \right.\\ \\[-7pt]{V_{(k,\,j-1)}} = {V_{(k,\,j - 1)}} + \left( {T_{(k,\,j - 1)}\cos {\alpha _{\left( {k,\,j - 1} \right)}} - {D_{\left( {k,\,j - 1} \right)}} - {m_{\left({k,\,j - 1} \right)}}g\cdot\sin {\gamma _{(k,\,j - 1,i,\,j)}}} \right)\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \\ \\[-7pt]{m_{(k,\,j - 1)}} = \left\{ {\begin{array}{*{20}{c}}{{m_{(k,\,j - 1)}} - {{\dot m}_1}\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}}}&{if}&{{t_{(k,\,j - 1)}}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_1}}\\\\[-7pt] {{m_{(k,\,j - 1)}} - {{\dot m}_2}\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}}}&{if}&{{t_1} < {t_{(k,\,j - 1) }}+\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}} \le {t_2}}\\ \\[-7pt]{{m_{(k,\,j - 1)}}}&{}&{otherwise}\end{array}} \right. \end{array}\end{equation}$

where $$ T_1 $$ and $$ T_2 $$ are the input thrust in the first and second stage, respectively, while $$ \dot m_1 $$ and $$ \dot m_2 $$ are the mass burning rate in the first and second stage, respectively, as computed using Equation (4). Note that the discrete non-dimensional objective functions using the discrete model formulation described above may be rewritten as

(25) $ \begin{equation}\begin{array}{l}C_{[(k,\,j - 1),(i,\,j)]}^{t} = \Delta {t_{(k,\,j - 1,i,\,j)}} = \frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}}\\\\C_{[(k,\,j - 1),(i,\,j)]}^{v} = {V_f} - {V_{(k,\,j - 1)}}\\\\C_{[(k,\,j - 1),(i,\,j)]}^{\gamma} = {\gamma _f} - {\gamma _{(k,\,j - 1)}}\end{array}\end{equation}$

Then, the total dimensionless cost functional for any arc $$ [(k,\,j-1),(i,\,j)] $$ which passes trough node ( $$ i,\,j $$ ) is

(26) $ \begin{equation}{\bar C_{[(k,\,j - 1),(i,\,j)]}} = {W_t}\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{\Delta {t_{\max }}.{V_{(k,\,j - 1)}}}} + {W_v}\frac{{{V_f} - {V_{(k,\,j - 1)}}}}{{\Delta {V_{\max }}}} + {W_\gamma}\frac{{{\gamma _f} - {\gamma _{(k,\,j - 1)}}}}{{\Delta {\gamma _{\max }}}},\end{equation}$

where

(27) $ \begin{equation}\begin{array}{l}\Delta {t_{\max }} = \mathop {max}\limits_k \,\frac{{\Delta {s_{(k,\,j - 1,i,\,j)}}}}{{{V_{(k,\,j - 1)}}}}\\[8pt] \Delta {V_{\max }} = \mathop {max}\limits_k \,|{V_f} - {V_{(k,\,j - 1)}}|\\[8pt] \Delta {\gamma _{\max }} = \mathop {max}\limits_k \,\,|{\gamma _f} - {\gamma _{(k,\,j - 1)}}|. \end{array}\end{equation}$

3.4 Integrating the problem formulation

The final CMCNF problem formulation is thus

(28) $ \fontsize{9.7}{10}\selectfont\begin{equation}\begin{array}{l}Minimise\,J = \sum\limits_{\{ (i,\,j) \in N\} } {\sum\limits_{\{ k \in K\} } {\left\{ {W_t}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{\Delta {t_{\max }}.{V_{(k,\,j \,{-}\, 1)}}}} \,{+}\, {W_v}\frac{{{V_f} \,{-}\, {V_{(k,\,j \,{-}\, 1)}}}}{{\Delta {V_{\max }}}} \,{+}\, {W_\gamma}\frac{{{\gamma _f} \,{-}\, {\gamma _{(k,\,j \,{-}\, 1)}}}}{{\Delta {\gamma _{\max }}}} \right\}{X_{k,\,j \,{-}\, 1,i,\,j}}} } \\ \\[-7pt] S.T.\,\left\{\! \begin{array}{l}\sum\limits_{\{ [(k,\,j \,{-}\, 1),(i,\,j)] \in A\} } {{X_{k,\,j \,{-}\, 1,i,\,j}}} \,{-}\, \sum\limits_{\{ [(k,\,j \,{-}\, 1),(i,\,j)] \in A\} } {{X_{i,\,j,k,\,j \,{-}\, 1}}} = \left\{{\begin{array}{*{20}{c}} {m \times n \,{-}\, 1} & \quad{for\,start\,node} \\ { \,{-}\, 1} & \quad {other\,nodes} \\\end{array}} \right. \\ \\[-7pt] {X_{k,\,j \,{-}\, 1,i,\,j}} \in \left\{ {0,1} \right\} \\ \\[-7pt] \Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}} = \sqrt {{{({r_{(i,\,j)}} \,{-}\, {r_{(k,\,j \,{-}\, 1)}})}^2} \,{+}\, {{({r_{(k,\,j\,{-}\,1)}}{\beta_{(k,\,j\,{-}\, 1)}})}^2}} \\ \\[-7pt] {\gamma _{(k,\,j \,{-}\, 1,i,\,j)}} = {\tan ^{ - 1}}\left( {\frac{{{r_{(i,\,j)}} \,{-}\, {r_{(k,\,j \,{-}\, 1)}}}}{{{{{{r_{(k,\,j\,{-}\,1)}} {\beta_{(k,\,j \,{-}\, 1)}}}}} }}} \right) \\ \\[-7pt] {{\dot \gamma }_{(k,\,j \,{-}\, 1)}} \cong \frac{{\Delta {\gamma _{(k,\,j \,{-}\, 1,i,\,j)}}}}{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{V_{(k,\,j \,{-}\, 1)}} \\ \\[-7pt] {L_{(k,\,j \,{-}\, 1)}} = {{m_{(k,\,j \,{-}\, 1)}}.{V_{(k,\,j \,{-}\, 1)}}.{{\dot \gamma }_{(k,\,j \,{-}\, 1)}} \,{+}\, {m_{(k,\,j \,{-}\, 1)}}.g.\cos {\gamma _{\left( {k,\,j \,{-}\, 1}\right) \to \left( {i,\,j} \right)}}} - {T_{(k,\,j \,{-}\, 1)}}\sin {\alpha _{(k,\,j \,{-}\, 1)trim}} \\ \\[-7pt] {\alpha _{(k,\,j \,{-}\, 1)trim}} = f(C{L_{(k,\,j \,{-}\, 1)}})\\ \\[-7pt] {CD _{(k,\,j \,{-}\, 1)trim}} = f(C{L_{(k,\,j \,{-}\, 1)}}) \\ \\[-7pt] \begin{array}{l}{T_{(k,\,j \,{-}\, 1)}} = \left\{ {\begin{array}{*{20}{c}}{{T_1}}&{if}&{{t_{(k,\,j \,{-}\, 1) \,{+}\, }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_1}}\\ \\[-7pt] {{T_2}}&{if}&{{t_1} < {t_{(k,\,j \,{-}\, 1) \,{+}\, }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_2}}\\ \\[-7pt] 0&{}&{otherwise}\end{array}} \right.\\\\[-7pt] {V_{(k,\,j\,{-}\,1)}} = {V_{(k,\,j \,{-}\, 1)}} \,{+}\, \left( {T_{(k,\,j \,{-}\, 1)}\cos {\alpha _{\left( {k,\,j \,{-}\, 1} \right)}} \,{-}\, {D_{\left( {k,\,j \,{-}\, 1} \right)}} \,{-}\, {m_{\left({k,\,j \,{-}\, 1} \right)}}g\cdot\sin {\gamma _{(k,\,j \,{-}\, 1,i,\,j)}}} \right)\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \\ \\[-7pt] {m_{(k,\,j \,{-}\, 1)}} = \left\{ {\begin{array}{*{20}{c}}{{m_{(k,\,j \,{-}\, 1)}} \,{-}\, {{\dot m}_1}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}}}&{if}&{{t_{(k,\,j \,{-}\, 1) + }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_1}}\\ \\[-7pt] {{m_{(k,\,j \,{-}\, 1)}} \,{-}\, {{\dot m}_2}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}}}&{if}&{{t_1} < {t_{(k,\,j \,{-}\, 1) + }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_2}}\\ \\[-7pt] {{m_{(k,\,j \,{-}\, 1)}}}&{}&{otherwise}\end{array}} \right. \end{array}\\ \alpha _{(k,\,j \,{-}\, 1)} = \left\{ {\begin{array}{*{20}{c}}{0}&{if}&{{t_{(k,\,j \,{-}\, 1) \,{+}\, }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_0}}\\ \\[-7pt] 4\alpha_{max} e^{-a[t_{(k,\,j \,{-}\, 1)}\,{-}\,t_0]}(e^{-a[t_{(k,\,j \,{-}\, 1)}\,{-}\,t_0]}\,{-}\,1)&{if}&{{t_0} < {t_{(k,\,j \,{-}\, 1) \,{+}\, }}\frac{{\Delta {s_{(k,\,j \,{-}\, 1,i,\,j)}}}}{{{V_{(k,\,j \,{-}\, 1)}}}} \le {t_1}}\\ \\[-7pt] 4\alpha_{max} e^{-a[t_{(k,\,j \,{-}\, 1)}\,{-}\,t_1]}(e^{-a[t_{(k,\,j \,{-}\, 1)}\,{-}\,t_1]}\,{-}\,1)&&otherwise\\\end{array}} \right.\\ \\[-7pt] \theta_{(k,\,j \,{-}\, 1)} =\alpha _{(k,\,j \,{-}\, 1)}\,{+}\,\gamma _{(k,\,j \,{-}\, 1,i,\,j)}\\ \\[-7pt] q\alpha_{(k,\,j \,{-}\, 1)} = 1/2\rho {V_{(k,\,j \,{-}\, 1)}^2}.{\alpha _{(k,\,j \,{-}\, 1,i,\,j)}}\\ \\[-7pt] \left| {q\alpha_{(k,\,j \,{-}\, 1)}} \right| \le q\alpha_{\max } \\ \end{array} \right. \\ \end{array}\end{equation}$

The right-hand side of the first constraint indicates that, if we want to determine optimum paths from the start node to every other node in the network, then in the minimum cost flow problem we set $$ b_{start} = m\times n - 1 $$ and $$ b_{i,\,j} = - 1 $$ for all other nodes. The minimum cost flow solution would then send unit flow from the start node to every other node along an optimum path.The flow parameter $$ X_{k,\,j - 1,i,\,j} $$ = 0 or 1 indicates whether each arc is in the feasible path ( $$ X_{k,\,j - 1,i,\,j}=1 $$ ) or not.

This formulation is used in the proposed algorithm described in the previous subsections, finally obtaining the optimum path. Since the flight parameters of this path are time dependent, the path is referred to as the trajectory.