This paper presents formulation and solution of long range flight vehicle and tactical air-to-air flight vehicle trajectory optimisation. The first case study is of a long range flight vehicle. Here an optimum steering program during powered phase has been evolved as control input for achieving maximum range with available propulsions in the presence of path and terminal constraints. The second case study is of a tactical flight vehicle for air-to-air application. Here a minimum flight time trajectory has been generated for covering a specified range pertaining to a specified air-to-air engagement by evolving pitch lateral acceleration as control input. Here also, there are many path and terminal constraints consisting of launch aircraft, pursuer, and evader. The studies have been carried out as part of system design activity of both flight vehicles. Both are real-life optimisation problems under several constraints. Through it is very difficult to solve such practical problems in flight dynamics using classical optimal control theory, it has been solved successfully using direct transcription method based on nonlinear programming. Rapid convergence has been achieved in four passes with minimum grids in first pass, to start with, and increasing the grids in subsequent passes. Solving such a real-life problem with proper convergence subjected to many constraints is claimed as novelty of present research.

NOMENCLATURE

 Ae Engine exhaust area CA,CN Flight vehicle axial, normal force coeffcient CD,CL Flight vehicle drag, lift force coeffcient D Drag J Performance Index L Lift M Mach Number, Mesh points Mmin(tb1+tc) Min Mach before second booster ignition (Nb,Nc) Number of grids during boost and coast Nb1,Nc1 Grid points during first booster burn, coast Nb2,Nc2 Grid points during second booster burn, coast P (.) Propagation operator Qdyn Dynamic Pressure $\left(\frac{1}{2}\text{ρ}{V}^{2}\right)\left(N/{m}^{2}\right)$ Rtogo Pursuer evader range to go S Pursuer aerodynamic reference area T Thrust Vm Pursuer velocity Vt Evader velocity (Vxa, Vya, Vza) Aircraft velocity components (east, north, up) (Vmx, Vmy,Vmz) Pursuer inertial velocity components (Vtx, Vty, Vtz) Evader inertial velocity components W Pursuer weight (mg) d Pursuer aerodynamic reference length (diameter) hk Discretization length at kth grid point mi ,m,mf Flight vehicle initial, instantaneous, final mass mp Propellant total mass ṁp Propellant burn rate ns Number of discretization segments (NLP) (ti, tf) Initial and final time (tc1, tc2) Coast time after first, second booster burnout u Control variable (x,y,z) Inertial position along (east, north, up) y State vector y(tf) State at final time Θ Parameters to be optimised in cost function (α, β) Angle of attack and side slip angle αr Resultant angle of attack η Load factor in ‘g’ (ηh, ην) Guidance demanded latax along yaw, pitch plane ω Actuator natural frequency (ϕ, γ) Azimuth, elevation angle along LV (ϕlma, γlma) Pursuer LOS wrt aircraft along yaw, pitch (ϕm, γm) Pursuer azimuth and elevation angle θ Flight vehicle attitude = (α + γ) ρ Air density ζ Actuator damping ratio

Subscripts

 a Mother aircraft h Horizontal (yaw) plane m Pursuer t Evader v Vertical (pitch) plane L Lower bound U Upper bound

Literature on trajectory optimisation of the flight vehicle (FV) is rather vast. The approaches on trajectory optimisation belong to two distinct categories1 such as (i) Direct method which uses mathematical programming by parameterizing the state and control histories and (ii) Indirect method which is based on solution of two-point boundary value problem (TPBVP) using optimal control principle. In essence, a direct method attempts to find the minimum of the objective function without constructing adjoint equations, control equations or transversality (boundary) conditions. An indirect method attempts to solve the optimal control state equations, boundary conditions, adjoint equations, control equations, transversality conditions and locates the roots of necessary conditions. The direct method appears to have been more popular in the application work over indirect method, possibly because of relative robustness1. Shrivastava & Reddy2 determined the optimal trajectory of a satellite launch vehicle under the constraint of product of dynamic pressure and angle of attack, control force and turning rate. Later Vathsal3, et al. solved the same problem using min-max approach. The fixed terminal time problem was solved using indirect method2,3. Hargraves4, et al. has used direct method to solve optimal evasion problem, namely, to obtain the optimal control of an aircraft being pursued by missile with proportional navigation guidance law such that the separation distance at the point of closest approach is maximised. Betts5 had carried out complete survey of numerical methods of trajectory optimisation using both direct and indirect methods. It is clear from literature survey that mathematical programming-based direct method, namely, NLP is more popular for solving real-life trajectory optimisation problems over indirect method.

First problem discussed here is long range FV range maximisation based on the literature survey. Ulybyshev6,7 developed a new technique to optimise continuous medium and low thrust orbit transfers. His approach combines large scale linear programming (LP) algorithms along with discretisation of trajectory dynamics on segments and a set of pseudoimpulses or control for each segment. Vavrina8, et al. developed a novel technique for low-thrust, interplanetary trajectory optimisation through hybridisation of a genetic algorithm and a gradient- based direct method to obtain optimal low-thrust trajectories for exploration of the solar system. Benson9, et al. carried out further research on trajectory optimisation using direct method where Gauss pseudo spectral method was used for orthogonal collocation of a NLP problem. Prasanna10, et al. carried out range maximisation of a real-life long range hypersonic research vehicle (HRV) using NLP. Second problem discussed is enhancement of launch-range (initial range between pursuer and evader) of a tactical FV air-to-air engagement. This was carried out by energy management of solid-propellant rocket motors. Menon11-12, et al. carried out synthesis of a high performance two-pulse-motor-propelled beyond visual range (BVR) vehicle by simultaneous optimisation of trajectory and propulsion system for air-to-air application using the direct method11,12. Kar13, et al. have evolved the maximum launch range pertaining to an air-to-air interception of a FV with one- pulse motor using NLP. The same methodology has been used to cover a given launch range of the present FV at minimum time subject to flight path and terminal constraints. Based on literature survey, it is evident that, trajectory optimisation of a practical problem in the presence of many constraints is handled by direct method. Present results on range maximisation of long range FV is based on the research carried out by Prasanna10, et al. In second problem, the minimum time trajectory for a given launch range pertaining to air-to- air interception of FV with two-pulse motor is extension of research carried out by Kar13, et al. which is also based on Prasanna10, et al. In both the cases, the optimisation has been carried out using NLP.

Based on literature survey it can be stated that since nineteen seventies, aerospace trajectory optimisation is still a fascinating field for applied research. Though a lot of open literature exists on space vehicle trajectory optimisation, the same for FVs under realistic path and terminal constraints is scanty. To fill this lacuna, the authors wish to alleviate it by present case studies. The first problem is of a long range FV. Here an optimum steering programme during powered phase has been evolved as control input for achieving maximum range with available propulsions in the presence of path and terminal constraints. The second problem consists of launch aircraft, pursuer, and evader for air-to-air interception. Here, pursuer being of BVR type, during midcourse guidance, the launch aircraft to track the evader. The formulation uses point mass kinematic model of launch aircraft, dynamical model of interceptor, and kinematic model of evader, and imposes several practical constraints in terms of flight path, terminal constraints, and line-of-sight (LOS) requirements. Evader information to be communicated by launch aircraft to the pursuer using data link during midcourse guidance. Also, data link maximum beam width is always constrained by on-board power requirement. So, pursuer has to be always within maximum beam width of launch aircraft (±30°) (Fig. 1) during midcourse guidance. Once the pursuer-evader range-to-go is less than 10 km, the terminal guidance starts using pursuer-mounted inboard seeker. During entire terminal guidance, the evader has to be always within the seeker maximum field of view for (±45°). These important practical constraints are based on hardware limitation and have to be modelled during synthesis of optimal pursuer trajectory. So the present study achieves realisation of optimal pursuer trajectory with inclusion of the following:

• Launch aircraft kinematic equations along with pursuer/ evader state dynamical equations in the formulation,
• Maximum line of sight (LOS) angle between pursuer and launch aircraft during midcourse guidance,
• Pursuer seeker gimbal freedom during terminal guidance as extra flight path constraints over and above other flight path and terminal constraints based on practical hardware considerations.

Figure 1. Air-to-air engagement (pursuer, evader and launch aircraft with γlma = 30).

2.1 Flight Vehicle Dynamical Model of Long Range FV

Governing dynamical equations of motion of FV assuming non-rotating spherical earth are (Fig. 3)

$\begin{array}{l}\stackrel{˙}{V}=\frac{T\mathrm{cos}\alpha -\frac{1}{2}\rho {V}^{2}S{C}_{D}}{m}-g\mathrm{sin}\gamma \hfill \\ \stackrel{˙}{\gamma }=\frac{T\mathrm{sin}\alpha +\frac{1}{2}\rho {V}^{2}S{C}_{L}\left(\alpha \right)}{m}-\left(\frac{g}{V}-\frac{V}{{R}_{e}+h}\right)\mathrm{cos}\gamma \hfill \\ \stackrel{˙}{X}=V\mathrm{cos}\gamma \hfill \\ \stackrel{˙}{Y}=V\mathrm{sin}\gamma \hfill \end{array}$(1)

$\text{here}\text{\hspace{0.17em}}\left(\alpha =\theta -\gamma \right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}g={g}_{0}{\left(\frac{{R}_{e}}{{R}_{e}+h}\right)}^{2}$(2)

where g0 = 9.81 m/s2 and earth radius Re = 6378137m. The initial conditions of state variables (V, γ, x, y) are known as (0 m/s, 90o, 0m, 0m) launch conditions for vertically launched FV. The control input is FV attitude θ time history. This has to be evolved through optimisation.

2.1.1 Cost Function

The PI used for maximisation of pursuer launch-range was

$\begin{array}{ll}J=\mathrm{max}{x}_{m}\left({t}_{f}\right)\hfill & \text{where}\text{\hspace{0.17em}}{t}_{f}={t}_{b}+{t}_{c}\hfill \end{array}$(3)

Here the parameter vector Θ = tc the element tc has to be maximised to obtain maximum range.

2.1.2 Flight Path Constraints

The constraints corresponding to flight path optimisation for a given propulsion system are as follow:

(i) The constraints on states, parameters, and control are

Figure 2. Free body diagram of present FV (pitch plane).

$\begin{array}{l}x\in \left(0,1500\right)\text{\hspace{0.17em}}\text{km}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}y\in \left(0,500\right)\right)\text{km}\hfill \\ V\in \left(0,4500\right)\text{\hspace{0.17em}}m/s\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\gamma \in \left(-90,90°\right)\hfill \\ {t}_{c}\in \left(0,500\right)s\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\theta \in \left(-90,90°\right)\hfill \end{array}$(4)

These numbers are completely based on probable maximum and minimum values of state variables in the problem.

(ii) Path constraint of FV being at attitude hold till 3.5 seconds after launch from the launcher.

$\begin{array}{l}\gamma \left(t\right)=90°\hfill \\ \left(0\le t\le 3.5\text{s}\right)\hfill \end{array}$(5)

(iii) Path constraint of FV after booster burnout is that during entire coast phase till impact the FV angle of attack should be zero.

$\begin{array}{cc}\theta -\gamma =0& \left({t}_{b}\le t\le {t}_{f}\text{s}\right)\end{array}$(6)

(iv) The terminal inequality constraints are

$\begin{array}{cc}M\left({t}_{f}\right)\ge {M}_{\mathrm{min}}\left({t}_{f}\right)& \text{where}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{\mathrm{min}}\end{array}\left({t}_{f}\right)=1.80$(7)

(v) Nonlinear state inequality constraints to be satisfied during complete trajectory as

$\begin{array}{l}-{\alpha }_{\mathrm{max}}\le \alpha \le {\alpha }_{\mathrm{max}}\hfill \\ 0\le \frac{1}{2}\rho {V}^{2}|\alpha |\le \frac{1}{2}\rho {V}^{2}{|\alpha |}_{\mathrm{max}}\hfill \\ -{\gamma }_{lma}^{\mathrm{max}}\le {\gamma }_{lma}\le {\gamma }_{lma}^{\mathrm{max}}\hfill \end{array}$(8)

here

In the present optimisation problem, α was considered as control input. The initial guess of the control history was arrived at based on burn out γ requirement for the class of FV considered in the present study14. Up to the power on phase, this control variable would be evolved based on specified PI [Eqn. (3)]. During free-flight phase the vehicle states were propagated till impact. The impact range was considered as the performance index to be maximised. The variable γ gets evolved through state equation [Eqn. (4)]. After convergence the attitude history is obtained as θ = α + γ. This optimized θ (t) till power on phase is a control input to the FV in real life.

2.2 Flight Dynamical Model of Tactical FV for Air- to-Air Engagement

2.2.1 Pursuer Dynamical Model

Governing pursuer dynamical equations of motion along pitch plane are

$\begin{array}{l}\frac{d{V}_{m}}{dt}=\frac{\left(T\mathrm{cos}\alpha -D\right)}{m}-g\mathrm{sin}{\gamma }_{m};\text{\hspace{0.17em}}\frac{d{\gamma }_{m}}{dt}=\left({\eta }_{vm}-\mathrm{cos}{\gamma }_{m}\right)\frac{g}{{V}_{m}};\hfill \\ \frac{d{h}_{m}}{dt}={V}_{m}\mathrm{sin}{\gamma }_{m}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\frac{d{x}_{m}}{dt}={V}_{m}\mathrm{cos}{\gamma }_{m}\hfill \end{array}$(9)

here ${\eta }_{m}=|{\eta }_{hm}|$. The initial conditions of state variables (Vm0, γm0, hm0, xm0 ) are available from engagement conditions.

The governing equations of motion of evader with constant velocity assumption are15

$\frac{d{V}_{t}}{dt}=0.0;\frac{d{\gamma }_{t}}{dt}={\eta }_{vt}\frac{d{h}_{t}}{dt}={V}_{t}\mathrm{sin}{\gamma }_{t}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\frac{d{x}_{t}}{dt}={V}_{t}\mathrm{cos}{\gamma }_{t}$(10)

The initial conditions of the state variables (Vt0, γt0, ht0, xt0) are available from engagement conditions.

2.2.3 Launch Aircraft Kinematic Model

The kinematic equations of motion of the launch aircraft are

${\stackrel{˙}{x}}_{a}={V}_{xa}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{\stackrel{˙}{h}}_{a}={\stackrel{˙}{V}}_{za}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\stackrel{˙}{V}}_{xa}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{\stackrel{˙}{V}}_{za}=0$(11)

Equation (11) shows that the aircraft after launch, will be at quasi-steady level flight with constant velocity. Initially pursuer and launch aircraft state variables are the same. At any instant of time, pursuer LOS wrt to aircraft can be written as (Fig. 2)

$\Delta x={x}_{m}-{x}_{a}\text{\hspace{0.17em}};\text{\hspace{0.17em}}\Delta h={h}_{m}-{h}_{a};\text{\hspace{0.17em}}{\gamma }_{lma}={\mathrm{tan}}^{-1}\frac{\Delta h}{\Delta x}$(12)

2.2.4 Seeker Gimbal Limit during Terminal Guidance

At any particular time, evader LOS angle wrt pursuer is

$\Delta {x}_{mt}={x}_{t}-{x}_{m};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\Delta {z}_{mt}={z}_{t}-{z}_{m};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda }_{p}={\mathrm{tan}}^{-1}\frac{\Delta {z}_{mt}}{\Delta {x}_{mt}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\theta }_{m}={\gamma }_{m}+a;\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda }_{p}-{\theta }_{m}$(13)

2.2.5 Cost Function

For minimising the total time of flight

$J=\mathrm{min}{\int }_{{t}_{i}}^{{t}_{f}\left(\Theta \right)}\text{dt}$(14)

This PI is used for pursuer to reach prespecified launch- range (viz., 110 km) at minimum time with available propulsion. The final time tf depends only on trajectory parameters.

So $\Theta ={\Theta }_{T}=\left({t}_{c1},{t}_{c2},{M}_{\mathrm{min}}\left({t}_{b1}+{t}_{c1}\right)\right)$(Fig. 3) and total parameter vector to be optimized consists of three parameters only.

2.2.6 Flight Path Constraints

The constraints corresponding to flight path optimisation for a given propulsion system are as follows:

(i) The constraints on states, parameters and control are

xm (0,150) km; hm (0,20) km; Vm ∈ (0,2500) m/s

Figure 3. Two-pulse motor firing sequence.

$\begin{array}{l}{x}_{m}\in \left(0,150\right)\text{\hspace{0.17em}}\text{km;}\text{\hspace{0.17em}}{h}_{m}\in \left(0,20\right)\text{\hspace{0.17em}}\text{km;}\text{\hspace{0.17em}}{V}_{m}\in \left(0,2500\right)\text{m/s}\hfill \\ {\gamma }_{m}\in \left(-90°,90°\right);\text{\hspace{0.17em}}{t}_{c1}\in \left(0,200\right)\text{s}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{t}_{c2}\in \left(0,200\right)\text{s}\hfill \end{array}$(15)

These numbers are completely based on interception geometry that is probable maximum and minimum values of state variables in the problem.

(ii) Path constraint of pursuer being at attitude hold till 2.5 seconds after launch from launch aircraft.

${\gamma }_{m}\left(t\right)=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(0\le t\le 2.5\text{s}\right)$(16)

(iii) The terminal inequality constraints are

$\begin{array}{ll}r\left({t}_{f}\right)\le {r}_{\mathrm{min}}\hfill & \hfill \text{where}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{\mathrm{min}}=5\text{m}\\ M\left({t}_{f}\right)\ge {M}_{\mathrm{min}}\left({t}_{f}\right)\hfill & \text{where}\text{\hspace{0.17em}}{M}_{\mathrm{min}}\left({t}_{f}\right)=0.80\hfill \end{array}$(17)

where

$r={r}_{{t}_{0}-{g}_{0}}=\sqrt{{\left({x}_{t}\left({t}_{f}\right)-{x}_{m}\left({t}_{f}\right)\right)}^{2}+{\left({z}_{t}\left({t}_{f}\right)-{z}_{m}\left({t}_{f}\right)\right)}^{2}}$

(iv) The nonlinear state inequality constraints are

(18)

$M\left({t}_{b1}+{t}_{c1}\right)\ge {M}_{\mathrm{min}}$ before 2nd booster start ${M}_{\mathrm{min}}\left({t}_{b1}+{t}_{c1}\right)=0.8$

The pursuer gimbal angle margin ${\lambda }_{p}-{\theta }_{m}={\lambda }_{gimbal}$, where ${\theta }_{m}=\alpha +{\gamma }_{m}$ and maximum gimbal angle limit is ${\lambda }_{gimbal}^{\mathrm{max}}=45°$.

3.1 Thrust, Mass and Aerodynamic Database

3.1.1 Thrust and Mass Data

Based on ground test of motor, vacuum thrust Tvac time history is available. At any particular time

$T\left(t\right)={T}_{vac}\left(t\right)-{p}_{h}{A}_{e}$(19)

Here ph is atmospheric pressure at given altitude. Ae is engine area with diameter of 0.85 m and 0.09 m for long range and tactical FV respectively. The long range FV initial mass mi is 12000 kg. At any instant $m={m}_{i}-{\stackrel{˙}{m}}_{p}t$. The final mass mf of the complete vehicle after power on phase is 3100 kg. After the first stage is jettisoned, the second stage mass for rest part of trajectory is 750 kg. The (mi, mf) of tactical FV are (165, 104) kg, respectively and maximum propulsion unit limit can be up to 61 kg.

3.1.2 Aerodynamic Database of Long Range FV
3.1.2.1 Aerodynamic Database of Complete FV

The wind tunnel database consisting of CA (M,h) is available. The normal force is calculated as

${C}_{N}={a}_{1}\alpha +{a}_{2}\alpha |\alpha |$ from ( a1 (M, h), a2 (M, h)) obtained from wind tunnel data base. Aerodynamic database needed for optimizing the PIs (Section 2) consists of CA (M, h) and CN (α, M) Cubic spline interpolation function csape from spline toolbox16 of MATLAB has been used to pre-process raw aero data ( CA, CN ). Spline toolbox has been used to fit a second degree polynomial over discrete zero database to ensure non zero value of second derivative to be used in optimisation routine for minimising PI using NR method for better convergence1. CL and C D from given CA and CN are

$\begin{array}{l}{C}_{L}={C}_{N}\mathrm{cos}\alpha -{C}_{A}\mathrm{sin}\alpha \hfill \\ {C}_{D}={C}_{N}\mathrm{sin}\alpha +{C}_{A}\mathrm{cos}\alpha \hfill \end{array}$(20)

Instantaneous angle of attack is evaluated as α = θ – γ. Based on (α, M, h) the aerodynamic coefficients ( CL , C D ) are computed from look up table database. The database of ( CL , C D ) has been obtained from ( CN, CA ) using Eqn. (20).

3.1.2.1 Aerodynamic Database of Second Stage of FV

The wind tunnel database consisting of a1 (M, h) and CA (M, h) are available for the second stage from a1 normal force is calculated as CN = a1α . Here also Spline toolbox has been used to fit a second degree polynomial over discrete aero database for better convergence. Calculation of CL and CD, from given CA and CN is similar to the complete vehicle as discussed. For all aerodynamic force calculations, reference area corresponding to 1 m diameter was used for both complete vehicle and second stage.

3.1.3 Aerodynamic Database of Tactical FV

In this case also, aerodynamic database needed for optimising the PI consists CA = (M, h) and CN (α, M). The cubic spline interpolation function csape of MATLAB has been used to carry out pre-processing of the raw aero data (CA , CN ). At any particular instant, ${C}_{L}={\eta }_{vm}/\left(\frac{1}{2}\text{ρ}{V}_{m}^{2}S\right)$ and corresponding (α, CD ) are computed from look up table obtained using Eqn. 20. For all aerodynamic force calculations, reference area corresponds 0.178 m diameter was used.

3.2 MATLAB Optimisation Toolbox

The function fmincon of optimisation toolbox is a general purpose performance index minimisation routine under linear / nonlinear equality and inequality constraints17. It finds the minimum of the problem specified by

(21)

where (x, b, beq, lb, ub) are vectors, (A, Aeq) are matrices, (c, ceq) are functions which return vectors and f (x) is a function that returns a scalar. The functions f (x), c(x) and ceq(x) can be nonlinear. While carrying out this exercise, the authors faced several practical issues which delineate this type of practical problem from the text book problems. Some of the practical issues corresponding to both the case studies are highlighted.

3.2.1 Case Study of Long Range FV

The variables used in optimisation were (x, y, V, γ, tcoast, θ). These were scaled to lie within [-1, 1] while using fmincon. Poor scaling can make a good algorithm bad. It changes the convergence rate, termination test, and numerical conditioning1. The optimisation was started initially in the first pass with 9 grid points during boost phase and 9 grid points during coast phase. The boost phase was for 51 seconds after which the coasting started. The coast phase was approximately for 500 seconds and tc had to be evolved through optimisation. In the present implementation total grids during boost Nb were chosen as odd numbers and total grids during coast Nc were the same as Nb . So in the first pass of trajectory optimisation there were total number of (9 + 9) = 18 grid points. In the second pass, Nb = 2 × 9 - 1 = 17 and Nc = 17 were chosen. The algorithm for choice of ( Nb, Nc ) for different iterations was and Nc = 17 were chosen. The algorithm for choice of ( Nb , Nc ) for different iterations is given. In the first pass, for total 18 grid points, the NLP constraints were formed using trapezoidal discretisation. In fourth pass, Nb = 65, Nc = 65 and total grid points were 130 and Hermite-Simpson discretization are used for discretisation. Total four passes were enough for achieving convergence in the present problem of interest. At this juncture, it is important to state that during boost phase, there was a rapid change in state and control variables especially (V, γ, θ). This calls for mesh refinement in different passes. The mess refinement algorithm is given.

3.2.2 Mesh Refinement Algorithm (Long Range FV)

First pass

Nb,1 = n where n is an odd number

Nc,1 = Nb,1

N1 = Nb,1 + Nc,1 = 2Nb,1

Second pass onwards

for i = 2,···, npass till convergence

Nb,i = 2Nb,i-1 –1

Nc,i = Nb,i

Ni = Nb,i + Nc,i = 2Nb,i

end

The optimisation was carried out only at the grid points and trajectories between the grid points was not known. To assess the accuracy of the results, after the convergence of optimisation process, the FV state dynamical equations were integrated using control histories obtained from interpolation of the NLP generated discrete controls θk. For NLP results to converge, tuning parameters such as maximum iterations, maximum function evaluations, tolerance for performance index convergence, factor for function derivative, and Hessian evaluation and tolerance for parameter convergence were chosen as (4000, 6000, 1.0e - 06, 1.0e - 06, 1.0e - 06).

3.2.3 Case Study of Tactical FV

While carrying out this exercise, the state variables were also scaled within [-1, 1] for better convergence of NLP. In the present problem, the pursuer trajectory with two-pulse motor configuration had four distinct zones, such as first boost, first coast, second boost, and second coast. The optimisation was started initially in the first pass with 5 grid points during the first zone, and 4 grid points during subsequent zones. In the present implementation, total grid points during first boost was chosen as odd number Nb1 and total grid points during other zones (Nc1, Nb2, Nc2) were chosen as even numbers (Nb1 - 1). So, in the first pass of trajectory optimisation there were total 17 grid points. In subsequently passes, mesh-refining was carried out according to following algorithm. Four passes were enough for achieving convergence in the present problem. In the fourth pass, (Nb1, Nc1, Nb2, Nc2) = (33, 32, 32, 32). NLP constraints were formed up to third pass using trapezoidal discretisation, and in fourth pass using Hermite-Simpson discretisation with ηh as control variables during entire mission. The discretization techniques are available in Betts1 and Arijit13 et al. and are not reported here for brevity.

3.2.4 Mesh Refinement Algorithm (Tactical FV)

First pass

Nb1,1 = n where n is an odd number

Nc1,1 = Nb1,1 –1; Nb2,1, = Nb1,1 –1 ; Nc2,1, = Nb1,1 –1

N1 = Nb1,1

Second pass onwards

for i = 2, …, npass till convergence

Nb1,i = 2Ni-1 –1

Nc1,i = Nb1,i –1; Nb2,i = Nb1,i –1; Nc2,i = Nb1,i –1;

Ni = Nb1,i

end

4.1 Case Study of Long Range FV

4.1.1 FV Optimal Trajectory Synthesis

In the present case study, for the given FV configuration along with the thrust, mass, and aerodynamic data, maximum range trajectory has been evolved under specified constraints. NLP converged at fourth pass with total grid points as 130. The control input θ time history till power on phase has been obtained as decision variable from optimisation at the grid points (Section 3). This has been used as input to point mass 2-degree of freedoms (2-DOFs) trajectory simulation model of the FV. Corresponding to the optimised trajectory, time history of FV downrange and altitude are shown in Fig. 4. Time history of (γ, θ, α) of the FV are shown in Fig. 5. Time history of FV Mach number is shown in Fig. 6. In all the figures, the kinematic variables obtained through point mass study match very well with the corresponding NLP outputs at grid points. This clearly demonstrates the correct convergence of optimisation algorithm along with input data compatibility.

Figure 4. FV downrange, altitude time history (130 grids +2-DOFs simulation with optimised θ till power on).

Figure 5. FV (γ,θ,α) (130 grids +2-DOFs simulation with optimised θ till power on).

Figure 6. FV Mach number (130 grids +2-DOFs simulation with optimised θ till power on phase).

4.1.2 Validation of Optimal Trajectory Using 6-DOFs Simulation Model

It is well known that the 2-DOFs plant model is well- suited for trajectory optimisation. In reality, the demanded θ profile, as evolved above has to be tracked by autopilot/ controller. Any practical FV controller works in the presence of finite system lag, actuator uncertainty, and limiters such as dead band, backlash, and maximum allowable $\left(\stackrel{˙}{\delta },\text{\hspace{0.17em}}\stackrel{¨}{\delta }\right)$. Also notch filter is added in control loop to mitigate the undesired effect of structural frequencies. Additional filters are also used to remove the noise of sensors such as accelerometers and rate gyro. Inclusion of all these elements adds delay in the system. As a result, the maximum range Rmax as evolved by optimisation of 2-DOFs model is never achieved in reality. The 6-DOFs simulation with the optimised θ history has been carried out. It has been found the range achieved in 6-DOFs simulation is 97 per cent of Rmax obtained from 2-DOFs model which is as expected. In 6-DOFs model, control input θ(t) is tracked by autopilot. This generates the control surface deflections as actuator command at 10 ms interval. Commanded fin deflections are passed through four independent actuators to achieve the required deflections. The actuators have been modelled as a second-order system with command input/output transfer function as $\frac{\left({\delta }_{0}\right)}{{\delta }_{1}}=\frac{{\omega }_{a}^{2}}{{s}^{2}+2{\zeta }_{a}{\omega }_{a}s+{\omega }_{a}^{2}}\text{with}\text{\hspace{0.17em}}{\zeta }_{a}=0.7,\text{\hspace{0.17em}}{\omega }_{a}=7\text{\hspace{0.17em}}\text{Hz}$. The actuator nonlinearities consist of dead zone and backlash of 0.23 deg and 0.115 deg half-width, respectively. Time history of FV acceleration components and body rate components during ascend phase till power on phase are shown in Fig. 7 and Fig. 8, respectively. All the acceleration and body rate components are well within the system limits. Time history of angle of attack and side slip angle till power on are shown in Fig. 9 which are within aerodynamic design limit. It is important to state that aerodynamic plays role only during ascend phase. The time history of dynamic pressure (Q), total angle of attack αT and Q × αT are shown in Fig. 10 and these parameters are all within structural design limits. So through realistic 6-DOFs simulation of FV dynamics it was infer that 97 per cent of Rmax is achievable which satisfies all design constraints.

Figure 7. FV (αx, αy, αz) during ascent phase (6-DOFs simulation with optimised θ till power on).

Figure 8. FV body rates (p,q,r) (6-DOFs simulation with optimised θ till power on).

Figure 9. FV (α,β) (6-DOFs simulation with optimised θ till power on).

Figure 10. FV (qdyn, αtot, qdyn X αtot) (6-DOFs simulation with optimised θ till power on).

4.2 Case Study of Tactical FV

A specific engagement condition chosen as design point corresponding to present FV is described. The pursuer, initially at 8 km altitude and at supersonic Mach number 1.2 had to intercept head on non-manouevering subsonic evader initially at 6 km altitude, flying at constant Mach number 0.7. Our aim was to evolve feasible minimum time optimal trajectory for the above engagement subjected to satisfying all constraints (Section 2) to achieve the specified launch range (110 Km). Due to two pulse motors in the FV, the initial and final longitudinal acceleration of the configuration is (10.5, 7.7) g. The burn durations (tb1, tb2) of pulse-1 and pulse-2 motors are of both configurations are (4.12, 6.90) seconds, respectively. The propulsion burn rate of both motors are $\left({\stackrel{˙}{m}}_{p1},\text{\hspace{0.17em}}{\stackrel{˙}{m}}_{p2}\right)$ are (7.7, 3.9) kg/s. The present configuration has been used for range extension with single-pulse rocket motor recently13. The present design study is with two-pulse rocket motor. The aim of the authors is to use the same configuration with both single-pulse and two-pulse motors for covering wide variety of launch ranges.

4.2.1 Pursuer Optimal Trajectory Synthesis

Minimum flight time trajectory for pre-specified launch- range has been evolved using present two-pulse motor. Corresponding trajectory parameters (tc1 , tc2, Mmin(tb1 + tc1)) are (6.9 s, 113.6 s, 1.6). Pursuer (h,γ,ηp) corresponding to present configuration as obtained through NLP convergence and through point mass simulation are shown in Fig. 11. From the Fig. 11 it is clear that optimal grid points evolved through NLP optimisation lie along the state trajectories obtained by solving dynamical equations with optimised parameters. This figure clearly indicates the proper convergence of NLP optimiser. The convergence has been achieved within four passes with mesh grid size of (33, 32, 32,32) at the fourth pass. Time history of pursuer γlma along with pitch plane pursuer- evader trajectory is shown in Fig. 12. During entire time, we have γlma < 30° which satisfies the data link constraint. Time history comparison of pursuer (V, M) for the given configuration is shown in Fig. 13. Pursuer interception Mach no is 0.84 and the terminal constraint in Eqn. 7 is satisfied. Time history comparison of pursuer latex and resultant angle of attack along pitch plane is shown in Fig. 14. (ηp, α) corresponding to present configuration is (0.2 g, 5°).

In this paper, at first trajectory optimisation of a long range FV for range maximisation with different paths and terminal constraints has been described. At first θ profile till power on as control input evolved through optimisation has been used to validate optimisation results using 2-DOFs point mass model of the FV. Later same control input has been used in 6-DOFs simulation model to validate the point mass results. Later the second problem has been addressed pertaining to air-to-air engagement of a tactical FV Here using a two-pulse motor, a minimum time optimal trajectory has been synthesised corresponding to a specified launch range. Both the studies have been carried out as a part of system design. In both the cases, the optimisation has been carried out using NLP in an iterative manner. At present, solution of the same optimisation problem is being attempted using global optimisers such as genetic algorithm (GA) and ant colony optimisation (ACO) technique.

Authors thanked Dr N Prabhakar, CC (SAM) for encouraging present research. Mrs C Regy Joseph of ASL for generating 6-DOFs simulation results. ASL and DRDL authorities also for granting permission to publish present research on system design activity.

1. Betts, J.T. Practical methods for optimal control using nonlinear programming. SIAM, 2001.

2. Shrivastava, S.K. & Reddy, N.M. Design of optimal trajectory under design constraints for a satellite launch vehicle. Acta Astronautica, 1976, 3(3), 333–347. doi: 10.1016/0094-5765(76)90140-5

3. Vathsal, S.; Menon, P.K. & Swaminathan, R. Minmax approach to trajectory optimisation of multistage rocket vehicles. IEEE Trans. Aero. Elect. Sys., AES, 1977, 13(2), 179–187. doi: 10.1109/TAES.1977.308454

4. Hargraves, C.R.; Johnson, F.; Paris, S.W. & Rettie, I. Numerical computation of optimal atmospheric trajectories. J. Guid., Contr. Dyn., 1981, 4(4), 406–414. doi: 10.2514/3.56093

5. Betts, J.T. Survey of numerical methods for trajectory optimisation. J. Guid., Contr. Dyn., 1998, 21(5), 193–207. doi: 10.2514/2.4231

6. Ulybyshev, Y. Continuous thrust orbit transfer optimisation using large-scale linear programming. J. Guid., Contr. Dyn., 2007, 30(2), 427–436. doi: 10.2514/1.22642

7. Ulybyshev, Y. Spacecraft trajectory optimisation based on discrete sets of pseudo-impulses. In Proceedings of AIAA Guidance, Navigation and Control Conference and Exhibit, Honolulu, Hawaii, August 2008, AIAA-2008- 6276-CP,. doi: 10.2514/6.2008-6276

8. Vavrina, M.A. & Howell, K.C. Global low-thrust trajectory optimisation through hybridisation of a genetic algorithm and a direct method. In Proceedings of AIAA/ AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, August 2008, AIAA-2008-6614-CP. doi: 10.2514/6.2008-6614

9. Benson, D. A.; Huntington, G.T.; Thorvaldsen, T.P. & Rao, A.V. Direct trajectory optimisation and costate estimation via an orthogonal collocation method. In AIAA Guidance, Navigation and Control Conference and Exibit, Keystone, Colorado, August, 2006. doi: 10.2514/6.2006-6358

10. Prasanna, H.M.; Ghose, D.; Bhat M.S.; Bhattacharya, C. & Umakant, J. Interpolation-aware trajectory optimisation for a hypersonic vehicle using nonlinear programming. In Proceedings of AIAA Guidance, Navigation and Control Conference and Exibit, San Francisco, California, August 2005, AIAA-2005-6063-CP. doi: 10.2514/6.2005-6063

11. Menon, P. K. & Lehman, L.L. A parallel quasi-linearisation for air vehicle trajectory optimisation. J. Guid., Contr. Dyn., 1986, 9(1), 119–121. doi: 10.2514/3.20076

12. Menon, PK; Cheng, V.H.L; Lin, C.A, & Briggs, M. M. High-performance missile synthesis with trajectory and propulsion system optimisation. J. Space. Rock., 1987, 24(6), 552–557. doi: 10.2514/3.25952

13. Kar, P.K.; Mukherjee, Arijit; Sarkar, A.K. & Padhi, Radhakant. Range extension of an air-to-air engagement by offine and online trajectory optimisation. AIAA Atmospheric Flight Mechanics Conference 0-11 August 2011, Porland, Oregon. doi: 10.2514/6.2011-6339

14. Siouris, G M. Missile guidance and control systems. Springer, Ed 1, 2004.

15. Zarchan, P. Tactical and strategic missile guidance. AIAA Inc., Ed 4, 2002.

16. Carl, de Boor. Spline ToolboxTM 3 User Guide. The Mathworks, Ed 11, 2008.

17. MATLAB. Optimisation ToolboxTM 4 User Guide. The Mathworks, Ed 8, 2008.

 Mr M. Manickavasagam did his graduation and post graduation from MIT (Anna University) in Aeronautical Engineering and IISc Bengaluru in Aerospace engineering. He has been working as scientist at DRDL, Hyderabad since 1991 and now at ASL. His current research interests are system design for long range flight vehicle, simulation and modelling, trajectory optimisation and multidisciplinary design optimisation for flight vehicle. Mr A.K. sarkar did his graduation and post graduation from IIT Kharagpur, in Aeronautical Engineering and IIT Madras, in Computer Science. He has been working as scientist at DRDL, Hyderabad since 1984. While at DRDL, he pursued Doctoral study in engineering from IISc Bengaluru. His current research interests are flight mechanics, simulation and modelling, state and parameter estimation, guidance and optimisation techniques for flight vehicle system design. Dr V. Vaithiyanathan is Associate Dean and Professor at the School of Computing from SASTRA University in Thanjavur. His research interests are soft computing, image processing and cryptography.