NOMENCLATURE
- ACO
ant Colony Optimisation
$$ A_e $$
exit area of the engine nozzle,
$$ m^2 $$
-
$$ C_{k,\,j-1,i,\,j} $$
the cost associated with arc
$$ [(k,\,j-1),(i,\,j)] $$
- CD
drag force coefficient
- CL
lift force coefficient
- D
drag force, N
- g
gravitational acceleration,
$$ m/s^2 $$
- GA
genetic Algorithm
- (i, j)
network nodes
- h
altitude, m
-
$$ I_{sp} $$
specific Impulse
- L
lift force, N
- m
mass of SLV, kg
- MCNF
minimum Cost Network Flow
- MILP
mixed Integer Linear Programming
- MOCO
multi-objective Constrained Optimisation
-
$$ n_z $$
vertical load factor
- PSO
particle Swarm Optimisation
- p(i, j)
the label of node (
$$ i,\,j $$
)- q
dynamic pressure,
$$ N/m^2 $$
- r
sLV distance from Earth centre, m
-
$$ R_0 $$
mean radius of Earth, m
- s
distance (length), m
- S
sLV reference area,
$$ m^2 $$
- SLV
space Launch Vehicle
- T
thrust force, N
- V
true Velocity, m/s
- x
vector of design variables
-
$$ X_{k,\,\,j-1,i,\,\,j} $$
amount of flow on arc
$$ [(k,\,j-1),(i,\,j)] $$
Greek symbols
-
$$ \alpha $$
angle-of-attack, deg
-
$$ \beta $$
range angle, deg
-
$$ \gamma $$
flight path angle, deg
-
$$ \theta $$
pitch angle, deg
-
$$ \rho $$
air density,
$$ kg/m^3 $$
1.0 INTRODUCTION
The problem of trajectory optimisation is important in all space missions. The solution to this problem enables one to specify the optimum trajectory program that should be followed to achieve a specified mission objective while simultaneously satisfying the constraints.
The main motivation behind this work is to develop an algorithm for the off-line trajectory design of SLVs in the shortest possible time. Since this is a nonlinear optimisation problem, the requirements of the problem have caused it to be placed in the category of hard problems. Indeed, it is NP-hard(Reference Hromkovic1).
Some researchers have formulated the constrained nonlinear optimisation problem as an optimal control problem in order either to achieve an analytical (closed-form) solution or to use numerical approaches such as indirect shooting methods(Reference Lu and Pan6,Reference Gath and Calise13) . However, this approach is generally problematic for the following reasons: (1) the adjoint (co-state) equations increase the size and complexity of the problem, (2) the convergence of the corresponding iterative solution procedure is sensitive to the accuracy of the initial guess for the input variables and (3) the optimiser may stop at a local minimum of the function that is not a global minimum. To overcome this limitation, many attempts with different initial guesses may be needed. Many authors have investigated alternative approaches to the solution of the optimal control problem(Reference Ricciardi, Vasile, Toso and Maddockx8,Reference Lu, Wang and Liu14,Reference Ricciardi and Vasile15) .
Heuristic (meta-heuristic) algorithms are widely used for constrained nonlinear optimisation: Dileep et al.(Reference Dileep, Kamath and Nair4), a Particle Swarm Optimisation (PSO) technique was applied to the multi-objective optimisation problem, and, in Bairstow et al.(Reference Bairstow, Weck and Sobieszczanski-Sobieski16), an attempt was made to use a self-adaptive Differential Evolution (EA) algorithm. In Buontempo(Reference Buontempo17) and Roshanian et al.(Reference Roshanian, Bataleblu and Ebrahimi18), the Genetic Algorithm (GA) is used to tackle the mentioned problem. In Karsli and Tekinalp(Reference Karsli and Tekinalp19), the optimisation is conducted using a Simulated Annealing (SA) algorithm. However, meta-heuristics do not guarantee that an optimal solution will be found, because they do not have an analytical proof(Reference Hromkovic1) but produce a near-optimal solution. Also, due to the random nature of the search space in these methods, a different solution is obtained by each execution of the algorithm. Meanwhile, in the algorithm proposed herein, an optimal and exact solution is obtained only once passing through the solution space.
In Zardashti et al.(Reference Zardashti and Bagherian20), a new model for the constrained optimal trajectory planning problem based on a layered network flow structure was presented, followed in Zardashti et al.(Reference Zardashti, Nikkhah and Yazdanpanah21) by a multi-objective optimisation-based layered network flow structure.
In this paper, a modified layered network flow structure that aims to construct a nonlinear optimal and constrained trajectory for an SLV mission is applied. By constructing a network grid around the Earth and utilising the discrete dynamic equations of motion, the trajectory of the LV from the launch site to a Low Earth Orbit (LEO) is obtained exactly by minimisation of a cost functional subject to the nonlinear dynamics and mission constraints of the SLV. For this purpose, the remainder of this manuscript is organised as follows: In Section 2, the flight mechanics formulation of the problem, with the cost functional and constraints, is presented. In Section 3, the discrete mathematical model and its network flow structure are constructed utilising a constrained MCNF algorithm based on the optimality conditions. Numerical simulations are presented in Section 4. Finally, we summarise and conclude the paper in Section 5.
2.0 FLIGHT MECHANICS FORMULATION
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)
$ \begin{eqnarray}m\dot V = T\cos\alpha-D-mg\sin\gamma \end{eqnarray}$
$ \begin{eqnarray}mV\dot \gamma =T\sin\alpha+L-mg\cos\gamma+mV^2/r\cos\gamma\end{eqnarray}$
$ \begin{eqnarray}\dot m = -\frac {{T}}{{I_{sp}.g_0}} \end{eqnarray}$
$ \begin{eqnarray} \dot r = V\sin \gamma \end{eqnarray}$
$ \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:
$ \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
$ \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:
$ \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:
$ \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
$ \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:
$ \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):
$ \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.0 ALGORITHM DEFINITIONS AND MODELLING
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):
$ \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.1 The feasibility conditions
In this application, for each node k in the 
$$ (\,j-1)^{th} $$
layer, 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 $$
), whereas otherwise 
$$ X_{k,\,j-1,i,\,j}=1 $$
. By this definition, the feasible space condition for each node in the network is provided.
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)] $$
,
$ \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
$ \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:
$ \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.2.2 Determining the optimum path
After computing the optimal cost label of the goal node, the optimal path can be obtained using a backtracking approachFootnote 1 as follows: For instance, let p(it, jt) be the optimal cost label of the goal node. If 
$$ (l,\,j\ t-1) $$
is its parent (predecessor) node 
$$ (l,\,jt - 1) = PARENT(it,\,jt) $$
, then the arc 
$$ \left[ {\left( {l,\,jt - 1} \right),\left({it,\,jt} \right)} \right] $$
is on the optimal path. Continuing this backward scheme, the start node is finally reached, then all the waypoints on the optimal path are extracted.
This form is appropriate for the application of the layered network algorithm, which is an efficient single-pass algorithm. This means that the proposed algorithm can find the optimal path by passing through the solution space once, as opposed to methods such as optimal control and meta-heuristic approaches, which may require multiple passes over the entire state space.
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:
$ \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
$ \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
$ \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
$ \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
$ \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
$ \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
$ \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:
$ \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
$ \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
$ \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
$ \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
$ \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.
4.0 SIMULATION RESULTS
To verify the performance of the optimal trajectory algorithm discussed above, numerical simulations are performed. The algorithmic engine is based on MATLAB. All the experiments are performed on a 2.50-GHz Intel Core-i5 CPU with 4 GB of physical RAM running 64-bit Microsoft Windows 10.
A typical Earth-centred polar coordinate system is used as input to the algorithm. Figure 1 shows the digitised heterogeneous model of this area loaded into MATLAB. The resolution of this map is 800,000 nodes (
$$ m=4,000 , n=200 $$
) with 
$$ \Delta r $$
varying linearly from 5m (near Earth) to 1000m (far from Earth), and 
$$ \Delta\beta $$
varying linearly from 
$$ 0.00001\deg $$
(near the launch point) to 
$$ 0.01\deg $$
(far from it) using
$ \begin{align}\Delta \beta_i &= \frac{\Delta \beta_2-\Delta \beta_1}{m}(i-1);\Delta \beta_2=0.01;\Delta \beta_1=0.00001 \nonumber\\ \nonumber\\[-9pt]
\beta_{(i,\,j)}&=\Delta \beta_i.(i-1).\frac{\pi}{180}\\ \nonumber\\[-9pt]
\Delta r_j&= \frac{\Delta r_2-\Delta r_1}{n}(\,j-1);\Delta r_2=2000;\Delta r_1=0.0\nonumber\\ \nonumber\\[-9pt]
r_{(i,\,j)}&=Re+\Delta r_j.(\,j-1)\nonumber\end{align}$
The SLV is launched vertically from the Earth with the initial conditions at 
$$ t=0 $$
s presented in Table 2. As can be seen, the initial velocity is calculated form 
$$ V_0=R_0\omega_E cos(\lambda_{Launch}) $$
, where 
$$ \omega_E=7.292115\times 10^{-5} rad/s $$
is the Earth rotational velocity and 
$$ \lambda_{Launch}=35 ^{\circ} $$
is the latitude of the launch site.
Table 2 Initial conditions
The final boundary conditions are specified in Table 3. As described in the table, the target orbit is a circular LEO at an altitude of 250km above mean sea level.
Table 3 Final conditions
The simulation was performed and the results compared with those obtained using the GA method. The GA can be considered as a representative evolutionary and meta-heuristic method. This method, which uses the same model with standard pitch programming, is developed with 50 populations and 100 generations.
Figures 2 and 3 show a comparison of the results of the two methods for the altitude and velocity versus flight time, respectively. It can be seen that both methods obtain similar trends, but with a difference of flight time of about 15s at the end of the trajectory. This means that the proposed method has a shorter total flight time, indicating better optimisation than the other method. Moreover, the execution time of the proposed algorithm is about 13min, while the execution time required for optimisation by the GA is about 1.4h. Therefore, from the processing time point of view, the proposed method is also preferable.
Figure 2. The time history of the flight height.
Figure 3. The time history of the flight velocity.
Figures 4–6 show the time histories of the SLV’s mass, angle-of-attack and angle-of-attack multiplied by dynamic pressure, respectively. It can be deduced that, when using both methods, the SLV mass changes are similar, but there are some differences in the other parameters. The angle-of-attack in the proposed model is found from the standard Equation (9), while in the GA model, it is obtained from differences of the pitch program data and the flight path angle. The dynamic pressure constraint is active for a small interval of time.
Figure 4. The time history of the flight mass of the SLV.
Figure 5. The time history of the angle-of-attack.
Figure 6. The time history of the angle-of-attack multiplied by the dynamic pressure.
The simulation results indicate that the final boundary conditions are fulfilled and the constrained parameters lie in their allowable ranges.
5.0 CONCLUSION
A layered network flow model is proposed to formulate the optimal trajectory design for an SLV. According to the optimality conditions of minimising the time of flight (or maximising the payload mass) simultaneously with orbital conditions, a variant of the constrained MCNF algorithm is developed to solve the problem regarding the dynamic constraints explicitly and the equations of motion of the SLV. The numerical simulation results show that, from the optimisation point of view, the results of the proposed method include a lower flight time, indicating better optimisation compared with the GA method. From the processing time point of view, the duration of the proposed algorithm is acceptable compared with the other method. At the same time, although the modelling and implementation of the proposed algorithm are slightly more complex than the GA method, it is worth using because of the inherent advantages mentioned above. Therefore, in future work, a Three-Dimensional (3D) network model of the SLV will be developed by considering the rotation of the Earth.
