Path generation for flying robots in mountainous regions

Abstract

In this paper, a new path planning method is proposed to resolve the problem of two-dimensional terrain following flight of flying robots in mountainous regions. The performance criteria considered for this mission design could include either the minimum vertical acceleration or the minimum flying time. To impose the terrain following/terrain avoidance constraints, various approaches such as least square method, Fourier series method, Gaussian estimation method, and Chebyshev orthogonal polynomial are explored. The resulting optimal control problem is discretized by employing a numerical technique namely direct collocation and then transformed into a nonlinear programming problem. The efficacy of the proposed method is demonstrated by extensive simulations, and particularly, it has been verified that this method is able to produce a solution that satisfies all hard constraints of the underlying problem.

Keywords

Terrain following terrain avoidance path planning direct collocation method optimal control

Introduction

The collision with ground obstacles and terrain near airport has always been one of the most critical issues in aviation accidents. Especially, the probability of this collision could be much higher in low-altitude flights, and therefore landing conditions have always been regarded as some of the most precarious flight conditions. A safe flight in low altitude is one of the stringent requirements in defining a mission for airplanes and flying robots. Furthermore, in low-altitude flights such as in mountainous regions, it is highly desirable to find a trajectory that guaranties the safety of flight and reaches the airplane to the destination by optimizing certain factors including time, fuel, vertical acceleration, and so forth. However, each of these factors is considered based on either the design specifications or the problem conditions. For instance, the optimization of the vertical acceleration in a flight path enables the reduction of the structural damages, which itself plays a pivotal role in the construction technology and cost. Moreover, the vertical acceleration factor impacts the passenger comfort in an airplane.¹

In recent years, flying robots are widely utilized for urban missions such as relief and rescue and traffic control. To navigate and control these flying vehicles, there are flight limitations associated with terrain following (TF) and terrain avoidance (TA). For instance, given that traffic control is essential for metropolises, one of the most preponderant challenges that each designer faces is high-rise buildings in the area. Since it is so common for flying robots to fly at low altitude, prompt decisions should be made because there is no ample time for deciding and generating a suitable flyable path. Therefore, the flight paths must be generated and optimized based on the mission, which requires a sufficient identification of the mission, flight vehicle dynamics, constraints, and flight obstacles.

In low-altitude missions, TF and TA maneuvers are both considered at the same time. TA is the automatic avoidance of airplanes from nearing any obstacles and TF is the adjustment of minimum altitude from the ground and the rapid passing of the trajectory in order to optimize the terrain following. In fact, low altitude and high velocity are the two main tactics of this maneuver. Also, TF and TA could be used either simultaneously or separately, depending on the mission and problem requirements. Generally, the scope of all issues pertaining to TF maneuver for flying robots includes^2,3:

Trajectory generation with TF purposes to enhance both survivability and mission effectiveness;

Control law development to follow the actual or model trajectory regarding the flight vehicle maneuvering capabilities;

Sensors blending for terrain data update;

Evaluation of the integrated system for different scenarios as well as emergencies.

The main goal in the trajectory design is obtaining an optimal path that strives to follow terrain in a way that all of the requirements, constraints, and expectations are met. Since this path is usually a curve, the TF problem may be simplified as a two-dimensional problem solved merely in a vertical plane perpendicular to the ground where the longitudinal dynamics of the flying robot is involved. Solving this problem is generally based on modeling the terrain with a suitable curve that is generated by gathering the terrain data and deploying mathematical methods for selecting the best curve.³

A myriad of studies with a variety of methods have been conducted on the subject of path planning for TF and TA maneuvers. Lu et al.^4,5 have used inverse dynamics⁴ and predictive control methods⁵ for optimal trajectory generation with the scope of TF for flight vehicles. In these papers, a point mass dynamic model has been utilized for designing a curve in a two-dimensional vertical plane with the purpose of optimizing the flight time and terrain. Note that the inverse dynamics method for generating the optimal path is used for offline purposes, whereas the predictive control method generally is used for online purposes.

Asseo⁶ used the steepest descent algorithm for path planning in low altitude in order to accomplish the TF/TA maneuvers. The main feature of this technique is to divide the three-dimensional problem into two subproblems associated with the horizontal and vertical planes. This simplification makes the algorithm usable for online purposes.⁶

Some studies^2,3,7 investigated the two- and three-dimensional optimal trajectory generation for TF/TA maneuvers of unmanned aerial vehicles. Malaek et al. utilized the equations of point mass motion for the dynamics of flying vehicle and Chebyshev orthogonal polynomial method for the terrain modeling in two dimensions. This paper also presents a three-dimensional path planning method for which the kinematics of flying vehicles is exerted in the equations of motion. One of the most outstanding features of that work is imposing comprehensive TF constraints such as external threats in motion equations of flying vehicle. It should be mentioned that since that method is categorized as an indirect method, it is effective in terms of accuracy and reliability but requires a great computational time.

Rahim and Malaek^8,9 perform optimal path planning for TF/TA maneuvers using fuzzy logic method for ground robots, submarines, and unmanned flying vehicles. In that work, to simulate the TF maneuver, terrain with various gradients has been used and the flight altitude and gradient of the terrain have been modeled with the fuzzy method. In Rahim and Malaek,⁹ the implementation of fuzzy method was done in two phases. The first phase includes offline navigation in which path planning is done through cost function, flight vehicle dynamics, and terrain data with the fuzzy method. The second phase includes online navigation in a way that flying vehicle follows the designed path by passing the terrain and is capable of altering the parameters of membership functions through fuzzy learning method and in conclusion, it implements three-dimensional TF/TA maneuvers simultaneously.

In Khademi et al.,¹⁰ a three-dimensional path planning problem for TF/TA purposes has been studied. In that research, the path planning problem is formulated as an optimal control problem, which is converted to a nonlinear programming problem (NLP) based on direct transcription methods from which a three-dimensional trajectory is obtained.

Many researches have been conducted in the area of optimal trajectory generation, which are categorized in ground, aerial, and marine applications.^11–14 To solve the resulting optimization problem, the existing techniques could be categorized into two types^15,16:

Indirect methods;

Direct methods.

The indirect methods are based on the calculus of variation, which converts the optimization trajectory problem into a Hamiltonian boundary value problem. The Hamiltonian boundary value problem needs to be solved via a numerical algorithm. In contrast, the direct methods cast the optimization trajectory problem as a NLP with bounded dimensions. In indirect methods, some auxiliary variables named costate variables are incorporated into the problem, and the state and control variables are calculated through those costate variables. In direct methods, the state and control variables are derived directly without the introduction of extra variables. Since most of the costate variables in indirect methods cannot be described physically, a satisfactory initialization of those variables may be cumbersome. Moreover, the equations are sometimes very sensitive to the variation of the costate variables and a small variation may lead to large deviations in the solution. These reasons point to the complexity of finding a suitable initial value and signify that direct methods might be more favorable than indirect methods for solving an optimization trajectory problem.¹⁷

Despite the differences in these two methods, the main idea behind them is to transform the optimization trajectory problem into a standard optimization problem (such as nonlinear programming). Unlike optimization trajectory problems, there is no explicit dynamics in standard optimization problems so the optimization variables are not time dependent. This trait allows the possibility of presentation in various solution methods for standard optimization problems. Each of these two main methods for solving the optimization trajectory problems is utilized based on its merits and demerits as well as designers’ approach. For instance, one common demerit in both methods is the long duration of the solution time and their inability for online usages. It is worth mentioning that there are many techniques for each category. For example, inverse dynamics and steepest descent approaches belong to indirect methods and direct collocation and pseudo-spectral approaches belong to direct methods.

In this paper, the direct collocation method as a particular type of the direct method is employed to solve the optimization trajectory generation problem for flying in earth boundary layer. Most of the literature on this problem focuses on indirect methods, which limit the flexibility of imposing hard constraints and lead to suboptimal solutions. For instance, a flight bound constraint is normally imposed through a soft penalty term in indirect methods, whereas this type of constraint can be easily handled by direct methods. By applying the direct collocation method to the TF problem, a NLP will be obtained. There is a broad set of algorithms for finding a local or near-optimal global solution of these problems, which may be found in numerous toolboxes including FMINCON, SNOPT, IPOPT, NPSOL, and SQPLAB. The reader could refer to Leyffer and Mahajan¹⁸ for a comparison of the features of these solvers.

The outline of the paper is as follows. First, the path planning problem is formulated as an optimal control problem in the next section, where TF and TA maneuvers are both taken into consideration. In “TF problem” section, the TF problem is described in great detail. In “Direct collocation method” section, the problem is discretized using direct collocation method. “Terrain modeling” section is dedicated to terrain modeling. Simulation results are presented in “Simulation results” section, and the conclusions of this research are drawn in the final section.

Problem formulation

As mentioned in the previous section, the objective of this work is to design an optimal flyable trajectory in such a way that certain nonlinear dynamic constraints and linear/nonlinear nondynamic constraints are satisfied and besides a cost function is optimized. Due to the main scope of the optimal control problem which is concentrated to find a history of control variables in such a way that minimize or maximize a performance index and satisfy the constraints of the system, controllers should be found which optimize the cost function and satisfy all of the constraints.

As already noted, the safe flyable motion in vicinity of the Earth in vertical plane is named TF, whereas the motion in horizontal plane is called TA. If the dynamics of the flying vehicle is incorporated in the problem, the complexity may become extremely high. The main cause of this issue is germane to physical constraints. Due to the complexity in the three-dimensional problem of TF/TA, we study the two-dimensional motion in vertical plane (TF) in this work. However, our proposed technique could be applied to the integrated TF/TA as well.

In this problem, it is supposed that the constraints consist of first-order ordinary differential equations representing the flying robot system dynamics, initial and final boundary conditions, and certain conditions relevant to state and control variables. In this section, the above-mentioned constraints and their formulations in the TF flight problem are presented. We could describe the system motion behavior as differential equations in a general form as follows¹⁹

\cdot ¯ x = f (x ¯, u ¯, t)

(1)

where

x ¯

is a

n \times 1

state vector,

u ¯

is a

m \times 1

control vector, t stands for time, and f is a given function. The dynamic equations of motion for a flying robot in a vertical plane with a flat ground model and under the assumption of a point mass movement are as follows^3,20

x \cdot = V \cos γ

(2)

y \cdot = V \sin γ

(3)

m \cdot = \frac{T}{I sp g}

(4)

V \cdot = \frac{T \cos α - D}{m} - g \sin γ

(5)

γ \cdot = \frac{T \sin α + L}{mV} - \frac{g}{V} \cos γ

(6)

where x and y are robot’s coordinates in vertical plane, V is the robot velocity, m is the flying robot mass, and γ is the flight path angle. The collection of all these variables defines the state vector as shown in Figure 1.

Figure 1.

Schematic of state and control variables for flying robot under the assumption of a point mass in vertical plane.

In equation (4), $I sp$ is specific impulse that indicates the engine power and reflects the intensity of the engine thrust and mass flow.²⁰ In equations (5) and (6), L and D are aerodynamic lift and drag forces, respectively, which are defined as follows

{L = qSC L D = qSC D

(7)

where q is the dynamic pressure, and S is the flying vehicle reference wing area which itself is a function of air density (ρ) and velocity. C_L and C_D are the lift and drag coefficients which are functions of angle of attach (α) and Much number (M). The propulsion force may be defined as follows

T = T max (M, y) \times EN (η)

(8)

where η denotes the engine power which is introduced as throttle and alters in the 0 ≤ η ≤ 1 range. Whenever η is 0, the engine is off, and whenever η is 1, the engine has the maximum thrust force. Furthermore, the symbol EN represents a nonlinear function of η. In equation (8),

T max

indicates the maximum thrust force for

η = 1

. Equations (2) to (6) are defined as the main flying motion equations, which describe the direction of the state variables, i.e. the direction of motion and arriving to next point in a small time interval dt based on their current position.

The main issue with these equations is that despite the initial conditions of state variables, they cannot be solved readily. The reason is that due to the presence of some unknown parameters in these equations, variations of these parameters play preponderant role in flight trajectory and state variables. These parameters are the control variables and their adjustment plays a fundamental role in the system performance. In this work, there are two control variables that are selected based upon the problem type and flying motion equations in a vertical plane.

The first control variable is the angle of attack, which is controlled through elevator control surfaces. The second control variable is the throttle setting, which is indicated by η. By considering the time history of the angle of attack and its derivatives during the time, the elevator’s history could be obtained. The control and state variables considered in this work are restricted to be bounded because of the flying robot’s physical and performance limitations, which in turn may imply that they could only vary in predefined ranges. These ranges may depend upon the dynamics of the flying robot, its desired performance, and also the environmental conditions. The following constraints are implemented for control variables

α min \leq α \leq α max

(9)

0 \leq η \leq 1

(10)

So far, we have introduced the flying robot dynamics constraints. Other constraints are related to the departure and terminal zones in the mission potential corridor, which could be classified as equality and inequality constraints.

In this work, the initial and terminal conditions are considered to be only equality constraints defined as follows¹⁷

Ψ ¯ (x ¯ 0, t 0) = 0

(11)

Ψ ¯ (x ¯ f, t f) = 0

(12)

where

t 0

is the initial time;

x ¯ 0

is the state vector at the initial time; t_f is the final time (that can be assumed to be fixed or free);

x ¯ f

is the state vector at the final time; and

Ψ ¯

is a vector consisting of p functions of the state variables, initial and final time.

In path planning problems, there are also some environmental constraints, such as obstacles, box constraints over the state and control variables. These constraints are classified into two categories of equality and inequality constraints, as noted below¹⁷:

Equality constraints

C ¯ eq (x ¯, u ¯, t) = 0

(13)

Inequality constraints

C ¯ ineq (x ¯, u ¯, t) \leq 0

(14)

where

C ¯ eq

is a

q \times 1

vector and

C ¯ ineq

is a

p \times 1

vector. External constraints in this work are altitude constraints in terrain form such that the flying robots keep a prespecified distance from the terrain and obstacles for collision avoidance purpose. These constraints are visualized in Figure 2 and imposed as follows

Terrain (x) + h min (m) \leq y (x)

(15)

Terrain (x) + h max (m) \geq y (x)

(16)

Figure 2.

Optimal trajectory with flight bound for TF maneuver.

The performance index is generally considered as follows¹⁹

J = Φ (x ¯ f, t f) + \int_{t 0}^{t f} (L (x ¯, u ¯, t) d t

(17)

where L and Φ are nonnegative functions. In a flying mission, the practical cost function is usually selected based upon the minimum possible time, the shortest desired trajectory, the TF, the minimum fuel consumption, and the minimum control effort or a combination of these performance indexes.

This study also may investigate designing some suitable cost functions which results in trajectories that are more convenient with physical characteristics and constraints of the flying vehicle mission.

The minimum-time cost function could be defined as follows

J = \int_{t 0}^{t f} d t

(18)

The shortest desired trajectory cost function is determined as

J = \int_{t 0}^{t f} \sqrt{y \cdot (t) 2 + x \cdot (t) 2} d t

(19)

The TF cost function is defined as

J = \int_{t 0}^{t f} {k (h (t) - h (t) terrain + (1 - k)) \times ((h (t) trim + h (t) clearance) - h (t))} d t

(20)

where k is a constant coefficient between 0 and 1,

h (t)

is the robot altitude from the reference surface,

h (t) terrain

is the terrain altitude, and

h (t) clearance

is the maximum allowable relative altitude of flying robot to the terrain surface. The minimum-fuel consumption cost function is defined as

J = - \int_{t 0}^{t f} m \cdot (t) d t

(21)

Another cost function that could be important for designers is the minimum vertical acceleration. In particular, it is desirable to reset the vertical acceleration in horizontal path to zero. As mentioned in “Introduction” section, this cost function is important for a number of reasons. First, by optimizing this cost function, the flying robot can withstand structural pressures, and second the convenience of passengers is assured in commercial airplanes. The above-mentioned cost function is one of the motivations behind this work, which could be defined as follows

J = \int_{t 0}^{t f} (y ·· (t)) 2 d t

(22)

where

y ··

is calculated according to equation (23) as follows

\frac{d 2 y}{dt 2} = \frac{d V}{dt} \sin γ + \frac{d γ}{dt} V \cos γ \frac{d 2 y}{dt 2} = (\frac{T \cos α - D}{m} - g \sin γ) \sin γ + ((\frac{T \sin α - L}{mV} - \frac{g}{V} \cos γ) V \cos γ)

(23)

TF problem

The main objective of the two-dimensional pure TF is to find a suitable trajectory that is parallel to the ground to satisfy the dynamics equations and optimize a performance index. Given that the trajectory should be parallel to the terrain and this may not be always possible, the TF problem is not pure anymore. Since it is difficult to follow the mathematical model of TF, a new method should be developed.

In these cases, generally two different methodologies could be utilized to solve the aforementioned problems. In the first approach, the computed control commands are modified in such a way that is suitable for the existing flying robot dynamics. The second approach attempts to create a virtual terrain profile based on the actual one, which is mathematically well behaved and is also compatible with the limitations of the flying robot. In this approach, we may define a suitable flight bound so that the flying robot is restricted to fly in this bound within a defined minimum and maximum altitude. This admissible flying corridor could be fine-tuned in various areas by the mathematical model of terrain, the dynamics of the flying vehicle, and the given performance index. Finally, it is the obligation of the designer to generate a reliable trajectory based on the predefined mission and aforementioned criteria. A pictorial sample of the flight bound is depicted in Figure 3 where, unlike Figure 2, this bound is not parallel to the ground and could alter based on the flying robot specifications and mission requirements.

Figure 3.

Flight corridor which could be defined based on the flying robot specifications and mission requirements.

This approach may be considered in cases that there is a trajectory with high steep and the flying vehicle dynamics is not capable of performing the pure TF mission due to possible structural damages and the designer should try a virtual terrain following strategy. In this approach, the flying vehicle may be forced to follow the virtual modified terrain instead of the actual terrain in an optimal manner.²

In what follows, we will find an optimal control strategy for TF purposes. The optimal trajectory will be generated by employing the direct collocation method and casting the problem as a NLP. In the next section, the direct collocation methods will be explained in great detail.

Direct collocation method

The trajectory of a flying robot can be accurately modeled as a collection of some discrete points in a continuous coordinate system from the initial time to the final time. More precisely, if the number of these identified points is sufficiently high, a suitable continuous trajectory can be recovered as a passable approximation of the main path. The approximate discrete and continuous trajectories have been shown in Figure 4.

Figure 4.

The approximate discrete and continuous trajectories.

The main idea behind direct collocation method is breaking up the state and control variables through a linear or higher order interpolation while satisfying all of the constraints. In this method, an optimal control problem is transformed into a NLP. Hence, this is an approximation technique that is employed as a numerical integration.^21,22

In this method, the state and control variables are divided into n sections and n+1 nodes as follows

t 0 = t 1 < t 2 < \dots < t i < \dots < t n + 1 = t f i = 1, 2, \dots, n + 1

(24)

where t₀ and t_f are the initial and terminal times, respectively. The state variable in node i, the control vector, the derivative of the state variable, and the cost value in each discrete interval between nodes i and i+1 are denoted as

x i ¯

u i ¯

f i ¯

, and J_i in Figure 5, respectively.

Figure 5.

Schematic of a discrete trajectory with state and control variables in each node and the cost function value in each discrete interval.

The time interval between nodes i and i+1 is equal to

h i = t i + 1 - t i

(25)

The approximation of the continuous trajectory by a set of nodes eliminates the differential equations of the flying system and reduces the problem to a nonlinear program. Many advantages of discretization methods have been discussed in the literature. Some of these discretization methods are Euler method, trapezoidal method, Hermite–Simpson method, and classic Runge–Kutta method.²² Two of these methods will be surveyed in this work in some detail.

Trapezoidal method

In this method, the difference between $x ¯ i$ and $x ¯ i + 1$ diminishes by increasing the number of nodes, leading to a smaller deviation from the true values. The error vector for each node is calculated as follows²²

x ¯ i + 1 - x ¯ i - \frac{h i}{2} [f ¯ i + f ¯ i + 1] = 0

(26)

Hermite-Simpson method

If Simpson integration rule and Hermite interpolation method are used instead of trapezoidal integration rule, the error vector for each node is calculated according to²²

x ¯ ci = \frac{1}{2} (x ¯ i + 1 - x ¯ i) - \frac{h i}{8} [f (x ¯ i, u ¯ i) - f (x ¯ i + 1, u ¯ i + 1)]

(27)

where

x ¯ ci

u ¯ ci

are the state and control variables vector in the centers of the ith and (i+1)th nodes, respectively, and are defined as

x ¯ ci = \frac{1}{2} (x ¯ i + 1 - x ¯ i) - \frac{h i}{8} [f (x ¯ i, u ¯ i) - f (x ¯ i + 1, u ¯ i + 1)]

(28)

u ¯ ci = \frac{1}{2} (u ¯ i + u ¯ i + 1)

(29)

In the above relations, i could change from 1 to n, and the vector d_i is referred to as the error in node i. If the state and control vectors in the nodal points are selected in such a way that the values of error vector are equal to 0 for $i = 1, 2, \dots,, n$ , the differential equations of motion are satisfied in a discrete form. Now, according to the discrete equations, which are a set of algebraic nonlinear equations, the NLP can be defined. It should be mentioned that this method results in the elimination of the differential equations, and therefore this method lays the foundations for the optimization process. A NLP could, in general, be presented with the variable $z \in R n$ such that the cost function

f (z)

(30)

is minimized subject to some algebraic constraints as follows^18,22

g (z) = 0 h (z) \leq 0

(31)

After the discretization of the equations through one of the aforementioned methods, the obtained NLP could be solved via a suitable numerical method.¹⁸

Terrain modeling

There are numerous methods to estimate the mathematical model of the terrain profile based on existing terrain discrete data. The least square method, Chebyshev orthogonal polynomial method, and Fourier series method are widely utilized for two-dimensional cases due to their good accuracy in estimation.²

The terrain modeling is performed in Figures 6 to 10 for a particular example using the aforementioned methods. The exhibited models in these figures are the results of the terrain mathematical model estimation by the least square method, Fourier series method, Gaussian estimation method, and Chebyshev orthogonal polynomial method, respectively. However, round-off error in these methods is more restrictive compared to Chebyshev orthogonal polynomial method. Increasing model accuracy while dealing with complicated terrains with irregular shapes and slopes is essential. In these estimations, we observed that the increase of the order of the approximation would be beneficial until a specific value is reached, and after that the increase would not result in improving the approximation and only causes the deflection of the curve in proportion to the true suitable model. Note that the determination of this suitable order is performed through various criteria, e.g. the error least square criterion.

Figure 6.

Terrain modeling through least square method.

Figure 7.

Terrain modeling through Fourier series method.

Figure 8.

Terrain modeling through Gaussian estimation method.

Figure 9.

Terrain modeling through Chebyshev orthogonal polynomial method.

Figure 10.

Comparison of the terrain models obtained from the above-mentioned methods.

Simulation results

Applying the constraints in TF problem

Optimal trajectory problems for a flying robot are generally defined over a period of time where the initial time is specified but terminal time could be specified or free based on the available constraints and performance index. In this case, we compute optimal trajectories influenced by the vehicle dynamics introduced as a constraint in the mathematical process.

The flying robot mathematical model for this research is taken from Khademi et al.¹⁰ Assume that the initial and final constraints of the state variables are specified as follows

X initial = {x (t i) = 0 (m), y (t i) = 1335 (m) M (t i) = 0.7, γ (t i) = 0 \circ, mass (t i) = 434.4 kg}

(32)

X final = {x (t f) = 7000 (m), y (t f) = 2600 (m) M (t f) = free, γ (t f) = free, mass (t f) = free}

(33)

The constraints originated from the flight path corridor are assumed to be modeled as the following inequality constraints for this case study

Terrain (x) + 30 (m) \leq y (x)

(34)

Terrain (x) + 350 (m) \geq y (x)

(35)

Also, a certain bound is considered for the state and control variables, which is imposed as inequality constraints with the specific values given in Table 1.

Table 1.

The physical constraints of the state variables for the flying robot.

Row	Physical constraints and their parameters	Value
1	Minimum angle of attack	$α min = - 3 \circ$
2	Maximum angle of attack	$α max = 12 \circ$
3	Minimum flight path angle	$γ min = - 80 \circ$
4	Maximum flight path angle	$γ max = 80 \circ$
5	The flying robot reference wing area	$S = 2.4999 m 2$
6	Primary mass of the flying robot	$m 0 = 434.4 kg$

In this case, since the final time of the flight mission is considered to be not specified but the initial and final positions of flying robot are specified, the system’s equations should be formulated in terms of x or y. The dynamics equations of the system are^2,3

\frac{dy}{dx} = \tan γ

(36)

\frac{d V}{dx} = \frac{T \cos α - D}{mV \cos γ} - \frac{g}{V} \tan γ

(37)

\frac{d γ}{dx} = \frac{T \sin α + L}{mV 2 \cos γ} - \frac{g}{V 2}

(38)

\frac{dm}{dx} = \frac{T}{I sp gV \cos γ}

(39)

Some physical features of the flying robot in this problem are indicated below

CL = α CL a

(40)

CD = CD 0 + KCL 2

(41)

CL a = 4.198 - (0.3425 M) + (1.0125 M 2)

(42)

Log (CD o) = - 2.02037 - (0.078043 Log (M)) + (0.05707 Log (y))

(43)

T = (η 6 + N (N, y)) T max (M, y)

(44)

where

CL a

and

CD 0

are provided in a tabulated form.³ As it is well known

CL α

is a function of Mach number

M

, while

CD 0

is a function of Mach number and flying altitude. In addition, the drag polar may have a more reliable form as it is usually needed for the purpose of flight test. It is interesting to note that, for the purpose of autopilot and control system design such a model may be advantageous as utilized in the proposed TF problem method.

Vertical acceleration minimization for a TF maneuver

As mentioned earlier, the main objective function considered in this work for designing a TF flight mode of the flying robot is concentrated to minimization of the vertical acceleration of the flight path with the aim of both avoiding some structural damages and supporting the main flight mission. In this problem, the initial and final conditions of the state variables are specified. The cost function for this problem is considered as follows

J = \int_{t 0}^{t f} (y ·· (t)) 2 dt

(45)

Given that the final time is unspecified in this example, the constraints are presented in terms of down range, x, and the cost function can be rewritten as follows

J = \int_{t 0}^{t f} (y ·· (t)) 2 d t = \int_{x 0}^{x f} ((\frac{T \sin α + L}{mV 2 \cos γ} - \frac{g}{V 2}) (1 + \tan 2 γ)) 2 d x

(46)

In fact, the goal is to find a suitable history of the control variables such as α, η through which the cost function is minimized.

Time minimization in vertical plane for TF

Apart from the vertical acceleration minimization problem for flying robot, the time optimization problem is another problem which is investigated in most researches. In this case, the cost function is minimized or maximized based on the problem’s requirements. Given the defined missions for flying robot, e.g. the urban traffic control, the cost function must be minimized and is defined as follows

J = \int_{t 0}^{t f} dt

(47)

Similar to the previous case, by changing the independent variable from t to x, the obtained cost function can be converted to equation (48)

J = \int_{t 0}^{t f} dt = \int_{x 0}^{x f} \frac{1}{V \cos γ} dx

(48)

The simulation results for both the minimum-time cost function and the minimum-vertical-acceleration cost function will be presented for two different methods of trapezoidal and Hermite–Simpson methods. For this purpose, the interval will be divided into 50 equal sections in these simulations.

The simulation results for the minimum-vertical-acceleration cost function are provided in Figures 11 to 16 including the variations of controls and states with respect to x.

Figure 11.

Throttle controller variations for optimal terrain following problem.

In Figures 11 and 12, variations of control variables including throttle and angle of attack versus downrange distance are shown.

Figure 12.

Angle of attack controller variations for optimal terrain following problem.

Figure 13 shows a typical case where the exact terrain could not be followed exactly due to the flying vehicle limitations. The flyable trajectory is therefore created in such a way that suits the flight dynamics and at the same time is minimized based on the cumulative vertical acceleration.

Figure 13.

Optimal trajectory with maximum altitude constraint.

In Figure 14, variations of flight path angle versus downrange distance are shown. In Figure 15, the flying robot Mach number history is presented. In Figure 16, variations of the flying vehicle mass versus downrange distance are shown.

Figure 14.

Optimal flight path angle.

Figure 15.

Optimal velocity in terms of Mach number.

Figure 16.

Mass variations of flying robot.

And finally, the flying vehicle actual trajectory that tracks the new path pattern instead of the exact terrain has been illustrated in the same figure.

In another case, we investigate a minimum time optimal trajectory during a TF flight. The simulation results for the minimum-time cost function using trapezoidal and Hermite–Simpson methods are given in Figures 17 to 22.

Figure 17.

Throttle controller variations for optimal terrain following problem.

Figure 18.

Angle of attack controller variations for optimal terrain following problem.

Figure 19.

Optimal trajectory with maximum altitude constraint.

Figure 20.

Optimal flight path angle.

Figure 21.

Optimal velocity in terms of Mach number.

Figure 22.

Mass variations of flying robot.

In this example, the flying vehicle specifications might not allow following the exact pattern of the terrain because of some structurally or aerodynamically restrictions. Therefore, it is necessary that the flying vehicle dynamics and structural characteristics be considered in the mathematical trajectory design approach.

In Figure 17, throttle setting histories related to this case are presented. As it is indicated later, the bang-bang scheme for throttle setting is based on the time optimal strategy. The bang-bang throttle setting history also indicates that the velocity, in general, will not be a constant during a TF flight. Switching times in this case depend on terrain profile and also the terminal conditions.

These simulations demonstrate that the obtained trajectory is near the terrain as much as possible and meets all of the constraints.

Conclusions and final discussions

The objective of this paper is to obtain an optimal reliable trajectory design approach to control a flying robot in the vicinity of the ground similar to TF flights. This approach concentrated on a two-dimensional vertical plane flight mode to satisfy the terrain constraints as well as the physical constraints of the flying robot. In particular, certain conditions on the dynamic variables, controllers, and initial and terminal boundary values are to be met.

This problem is formulated as an optimal control problem and then the obtained optimal control problem is transformed into a NLP through direct collocation method. The objective function of the optimization is also considered as the minimum time or the minimum vertical acceleration.

For the problem of designing an optimal trajectory for TF modes, one main challenge is the terrain modeling in a two-dimensional vertical plane. There are numerous methods for terrain modeling, four of which are utilized for terrain modeling in this work, including the least square method, Fourier series method, Gaussian estimation method, and Chebyshev orthogonal polynomial method. As observed in the paper, each one of these methods yields some errors with respect to the true model. Therefore, error least square criterion is utilized to select the best model and it is verified that the Fourier series method may have the smallest error in comparison to other methods, and as a result this method is used for the modeling in this paper. It should be noted that the challenging issue pertaining to the terrain modeling is that the created curve may pass under the terrain points in some positions, leading to a lower altitude than the terrain. Hence, the altitude should be compensated for through a coefficient to guarantee that the obtained trajectory will not collide with the terrain. A suitable approach for tackling this problem is to first calculate the vertical difference between the curve model and the true points, and then shift the curve upward in proportion to the maximum height difference through which a reliable trajectory can be generated.

After carrying out the terrain, a practical cost function is to be selected for meeting the required missions. The cost function considered in this paper is the minimum vertical acceleration, where a secondary goal is to reset the vertical acceleration in horizontal trajectory to zero. In other words, by optimizing the cost function, the flying robot undergoes less structural damages and the passenger comfort is assured in aerial transport vehicles. Moreover, since the minimum-time cost function is considered in some researches, this function is also investigated in this paper in addition to the minimum-vertical-acceleration cost function.

After formulating the path planning problem as an optimal control problem, direct collocation method is deployed. One of the most outstanding features of this method is that the differential equations can be simplified as algebraic equations. In this work, it is verified in simulations that after employing this method and solving the NLP, a feasible trajectory meeting all requirements of the problem and accomplishing the mission can be found. In other words, although the dynamic equations of the robot are highly nonlinear, the simulation results are very promising.

To utilize the direct collocation technique, trapezoidal and Hermite–Simpson methods are used in this work. The Hermite–Simpson method is a three-order approximation, whereas the trapezoidal method is a linear approximation. Although a higher approximation degree results in a more accurate solution, the computation time increases. It is observed that Hermite method is more accurate than the trapezoidal method, yet the computational complexity is higher. Moreover, these aforementioned methods can be chosen based on their merits, demerits, designers’ approach, defined missions, and so forth.

It should be mentioned that there is a great number of heuristic methods for solving the optimal path planning problem and most of them are not capable of solving this complex problem under tight constraints. The efficacy of direct collocation methods is elucidated in this work and it is demonstrated that the problem can be solved with a high accuracy. Therefore, this approach could be utilized in various optimization problems and this work could lay the foundations for further research in the future.

The development of a control system for special maneuvers is a quite lengthy process that demands careful attention in both software and hardware development in addition to numerous simulations and flight testing.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Appendix

References

Malaek

Abbasi

. Near-optimal terrain collision avoidance trajectories using elevation maps. IEEE Trans Aerosp Electron Syst 2011; 47: 2490–2501.

Malaek

Kosari

. Dynamic based cost functions for TF/TA flights. IEEE Trans Aerosp Electron Syst 2012; 48: 44–63.

Malaek SM and Kosari AR. A novel minimum time trajectory planning in terrain following flights. In: Aerospace Conference, Iran, 2003. Proceedings. 2003 IEEE, vol. 8, pp.8_3755–8_3762. Montana, USA: IEEE, 2003.

Pierson

. Optimal aircraft terrain-following analysis and trajectory generation. J Guidance Control Dyn 1995; 18: 555–560.

Pierson

. Aircraft terrain following based on a nonlinear continuous predictive control approach. J Guidance Control Dyn 1995; 18: 817–823.

Asseo SJ. Terrain following/terrain avoidance path optimization using the method of steepest descent. In: Aerospace and electronics conference, 1988. NAECON 1988, proceedings of the IEEE 1988 National, pp.1128–1136. Dayton OHIO: IEEE, 1988.

Malaek SM, Kosari AR and Jokar S. Dynamic based cost functions for TF/TA flights. In: Aerospace conference, 2005 IEEE, pp.1869–1876. Montana, USA: IEEE, 2005.

Rahim M and Malaek SM. Intelligent operation using terrain following flight in unmanned aerial vehicles. In: Aerospace conference, Big Sky, MT, 3 March 2007, pp.1–8.

Rahim

Malaek

. Aircraft terrain following flights based on fuzzy logic. Aircraft Eng Aerosp Technol 2011; 83: 94–104.

10.

Khademi I, Maleki B, and Mood AN. Optimal three dimensional terrain following/terrain avoidance for aircraft using direct transcription method. In: 2011 19th Mediterranean Conference on Control & Automation (MED), 20–23 June 2011, pp.254–258. Corfu, Greece: IEEE, 2011.

11.

Betts

. Survey of numerical methods for trajectory optimization. J Guidance Control Dyn 1998; 21: 193–207.

12.

Rao

. A survey of numerical methods for optimal control. Adv Astronaut Sci 2009; 135: 497–528.

13.

Conway

. A survey of methods available for the numerical optimization of continuous dynamic systems. J Optim Theory Appl 2012; 152: 271–306.

14.

Huang

Nan

. A survey of numerical algorithms for trajectory optimization of flight vehicles. Sci China Technol Sci 2012; 55: 2538–2560.

15.

Von Stryk

Bulirsch

. Direct and indirect methods for trajectory optimization. Ann Oper Res 1992; 37: 357–373.

16.

Grimm

Hiltmann

. Direct and indirect approach for real-time optimization of flight paths. Optimal control, Berlin Heidelberg: Springer, 1987, pp. 190–206.

17.

Dai R. Three-dimensional trajectory optimization in constrained airspace. PhD dissertation, Faculty of Auburn University, Alabama, 2007.

18.

Leyffer S and Mahajan A. Nonlinear constrained optimization: methods and software. Argonne, IL: Argonne National Laboratory.

19.

Kirk DE. Introduction to Optimal Control Theory. Englewood Cliffs, NJ: Prentice Hall, 1970.

20.

Tian B, Zong Q and Dou L. Ascent phase trajectory optimization for aircraft based on nonlinear programming. In: Control and Decision Conference (CCDC), 2010 Chinese, pp.2076–2081. Xuzhou, China: IEEE.

21.

Hargraves

Paris

. Direct trajectory optimization using nonlinear programming and collocation. J Guidance Control Dyn 1987; 10: 338–342.

22.

Betts JT. Practical methods for optimal control using nonlinear programming, volume 3 of Advances in Design and Control. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), 2001.