Trajectory optimization for parafoil delivery system considering complicated dynamic constraints in high-order model

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.1 (1-13)

Aerospace Science and Technology ••• (••••) ••••••

1

67

Contents lists available at ScienceDirect

68

2 3

69

Aerospace Science and Technology

4

70 71

5

72

6

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13

Trajectory optimization for parafoil delivery system considering complicated dynamic constraints in high-order model ✩

77 78 79 80

14 15 16 17 18 19

a

b

Hao Sun , Shuzhen Luo , Qinglin Sun

a,∗

a

a

, Zengqiang Chen , Wannan Wu , Jin Tao

c

a

22 23 24 25 26 27 28 29 30 31 32 33 34

82 83

College of Computer and Control Engineering, Nankai University, Tianjin, China b Mechanical and Aerospace Engineering, Rutgers University, NJ, United States c Aalto University, Helsinki, Finland

84 85 86

20 21

81

a r t i c l e

i n f o

Article history: Received 26 June 2019 Received in revised form 27 September 2019 Accepted 6 December 2019 Available online xxxx Keywords: Parafoil delivery system Trajectory optimization Gauss pseudospectral method Homing theory Dynamic constraint Multiple control objective

35 36 37

87

a b s t r a c t

88

Trajectory optimization is essential for all unmanned aerial vehicles, especially the parafoil delivery system, as the optimized trajectory must be practical and realizable. However, with the traditional 3-degree of freedom (DOF) model used in trajectory optimization, the realizability of the optimized trajectory is largely reduced without consideration of its dynamic constraints. Moreover, the landing and tracking precision of the system will also be reduced. To solve this problem, in this paper, a novel trajectory optimization method is explored. Different from existing research, trajectory optimization is achieved with a complete 6-DOF dynamic model of the parafoil delivery system. When this is combined with the Gauss pseudospectral method, complicated dynamic constraints can be taken into account in the trajectory optimization. In addition to achieving a global optimal control objective under dynamic constraints, the proposed method also can include multiple external environment constraints, such as terrain, wind disturbance and flared landing. Finally, our results show that the proposed method can achieve all the control objectives successfully. The optimized trajectory is smooth and realizable. Compared with other trajectory optimization methods that use a 3-DOF model, the advantage of the proposed method is clear. Considering the complicated dynamic constrains of the system, it has much better realizability and is more consistent with reality. © 2019 Elsevier Masson SAS. All rights reserved.

89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

38

104

39

105

40

106 107

41 42

1. Introduction

43 44 45 46 47 48 49 50 51 52 53 54 55 56

A parafoil delivery system is a special kind of unmanned aerial vehicle (UAV). This system consists of a parafoil canopy and a payload. The parafoil canopy will be full of air and provide the lift force during the flight process. Furthermore, the control system and the supplies are all set on the payload. The only horizontal control quantity of the system is the flap deflection, which is achieved by the deflection of the trailing edge of the parafoil canopy. The yaw angle of the system can be controlled by the flap deflection. This system is widely used in airdrop supply, aircraft recovery and other applications [1]. For example, a 4200 ft2 parafoil delivery system was applied as the recovery system for X-38 [2]. For all kinds of UAVs, the trajectory optimization is important, especially for the parafoil delivery system [3]. Regarding the appli-

57 58 59 60 61 62 63 64 65 66

✩

This work was supported in part by the National Natural Science Foundation of China under Grant (61973172) and Grant (61973175), and in part by the Key Technologies Research and Development Program of Tianjin under Grant (19JCZDJC32800). Corresponding author. E-mail address: [email protected] (Q. Sun).

*

https://doi.org/10.1016/j.ast.2019.105631 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

cation of the parafoil delivery system, homing control is the most important mission. The payload must be delivered to the target position accurately. Compared with accurate human-in-the-loop control [4], in the automatic homing process, trajectory optimization needs to be achieved before trajectory tracking. The optimized trajectory must be accurate and realizable. However, from studies of its aerodynamic force [5,6], considering the flexible shape of the parafoil canopy, the parafoil delivery system is a nonlinear system and has special dynamic characteristics. All of these dynamic characteristics form the dynamic constraints of the parafoil system; for example, with a nonlinear relationship between the flap deflection and the angular velocity of the yaw angle, this system has large inertia and cannot turn sharply; the flight velocity and angular velocity will obviously be changed by the wind disturbance [7,8]; not only the horizontal states, but also the vertical velocity and the accelerated velocity will be changed by the horizontal control quantity. These dynamic constraints place higher requirements on the trajectory optimization of the parafoil delivery system. In complex external environments, the nonlinear characteristics will be more distinct. Therefore, it is necessary to consider dynamic constraints and multiple environmental constraints in the trajectory optimization. Otherwise, the optimized trajectory may be unrealis-

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.2 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

tic, such as a control objective unreachable by flap deflection, an overlarge turning in a windy environment or being blocked by terrain. In recent years, there have been two main ideas involved in the studies on trajectory optimization: the traditional method and the optimal method. First, in both methods, the 3-DOF or 4-DOF model is most commonly applied for trajectory optimization [9–11]. The 3-DOF model has almost no consideration for the dynamic constraints of the system. It is mostly applied because of its faster computing speed, especially for in-flight optimization missions [12]. Previous studies have mostly employed high-order models, such as 6-DOF or 8-DOF models, to track the optimized trajectory [13,14]. In the 3-DOF model, both the horizontal and vertical flight velocities are constant. However, it is obvious that the flight velocity is time-varying with the flap deflection of the parafoil system and the external wind disturbance in both high-order models and actual experiments [15,16]. Therefore, there definitely are differences between the optimized trajectory and the actual trajectory. It is therefore necessary to consider the dynamic constraints of the parafoil delivery system in the trajectory optimization. However, in the optimization, the first challenge is to translate the highorder model into the form of a partial differential equation. In this way, the system states X = [x1 , x2 , . . . , xn ] T can be updated by X˙ = [˙x1 , x˙ 2 , . . . , x˙ n ] T . The main problem is to combine the nonlinear model and the optimization method by resolving the complex relationships among these states in the model. Furthermore, the high-order models are obviously nonlinear and much more complicated than the 3-DOF model. Therefore, the control quantity may have large oscillation in the optimization making the optimized trajectory difficult to track. In recent years, with the continued development of optimization methods, it is possible to achieve trajectory optimization with the high-order models by choosing a correct optimization method and an effective design of the objective function. This will be highly advantageous for special missions that have strict requirements regarding landing precision. By analyzing the existing research concerning these two ideas, it can be observed that the traditional method and the optimal method have different advantages and applications. The traditional method is designed specifically for actual airdrop experiments: For instance, a classic approach incorporating the traditional method was applied for the recovery system of X-38 [2,17,18]. This method has several phases, including the deployment, acquisition, homing, energy management, final stage, and landing stage. In this method, there are several types of optimized trajectories, such as Dubins path construction and T path construction. Studies on Dubins path include those by Rosich [13], Rademacher [19] and Yomchinda [20] and studies on the T path include those by Jann [3], Cleminson [21] and Culpepper [22]. Ref. [23] applied the Dubins path for homing control of the parafoil delivery system with a wind measurement. The results demonstrate reductions in the landing position error by more than 30% as well as a complete elimination of the possibility of downwind landings. Moreover, although there are different types of trajectories, all the optimized trajectories are composed of straight lines and standard circles. With the simple and standard trajectories, the optimized trajectory can be tracked in the airdrop experiment [24,25]. Moreover, the computing power requirement for the optimization algorithm is not very high for the 3-DOF model and these standard trajectories, which is also an important advantage. The optimization algorithm includes the CCRRT method from Luders [26], the band-limited guidance method from Carter [27], and other methods. The traditional method is widely applied in experiments because it is convenient for trajectory tracking. However, it obviously cannot consider the dynamic constraints of the system with the 3-DOF model. The actual flight path also cannot be an ideal straight line or standard circles in a windy environment. Our previous research [24] suggests that, in

a windy environment, the error between the reference trajectory and the actual flight trajectory is greater than 15 m. Ref. [21] also bluntly points out that the traditional method can render solutions more tractable and computationally efficient, albeit at the expense of some fidelity. The other approach to trajectory optimization is the optimal method. In this method, the optimized trajectory is an irregular curve. In the research area of parafoil delivery system trajectory optimization, a description of this method was presented by Gimadieva [28]. In this reference, the boundary-value problem is solved by the method of Krylov and Chernous’ko. A 3-DOF model has also been applied for trajectory optimization [29], such as a direct multiple shooting method by Babu [30] and a CC-RRT method by Luders [31]. Using this method, the irregular trajectory is difficult to directly track in actual experiments. Therefore, for tracking, the curve is often divided into several straight paths and is mostly tracked by a high-order model or a real system [14,29]. For the optimal method, the smoothness of the trajectory is very important. Moreover, with the development of computer technology, the optimal trajectory method has also been improved. In studies by Slegers and Rogers, a simplified 6-DOF model was applied to Monte Carlo simulation experiments concerning trajectory optimization [32,33]. The optimized trajectory was achieved with many repeat computations. Then, this method and another guidance method based on a graphics processing unit (GPU) were applied to flight experiments [34]. However, the computation time for this method with the Monte Carlo simulation experiments is very long. Furthermore, Ref. [32] points out that the total miss distance is partly due to poor geometry at the initiation of terminal guidance and partly due to unanticipated deadband in the innerloop of the horizontal controller. It is obvious that these errors are caused by ignoring dynamic constraints in the simplified model. The poor geometry at initiation of terminal guidance ignores the constraint of the angular velocity of the yaw angle, and the unanticipated deadband of the horizontal controller denotes neglect of the control constraint. All of these studies illustrate that, to achieve an accurate and realizable optimized trajectory, it is important to further consider the dynamic constraints of the parafoil delivery system. In this paper, a 6-DOF model of the parafoil delivery system is presented first. By considering the flight velocity and the angular velocity of the system, its dynamic constraints can be taken into account in the model. Meanwhile, a 3-DOF model is also outlined for comparison. Then, to improve the realizability of the optimized trajectory, a novel trajectory optimization method is designed based on the Gauss pseudospectral method (GPM). In this section, the theory of the Gauss pseudospectral method is explained. Twelve states are involved and updated in the optimization, including the position, velocity, Euler angle and angular velocity. In this way, utilizing the principle of the Gauss pseudospectral method, the complete 6-DOF model can be involved in the trajectory optimization. Not only the dynamic constraints of the flight velocity and angular velocity, but also the initial (terminal) constraint, control constraint and real-time path constraint are all considered. All the other constraints and the control objectives are introduced into the optimization, including the wind disturbance, terrain avoidance, and landing precision. Finally, detailed simulations are carried out with various initial positions and external environments. The results show that all the control objectives are achieved successfully and the landing precision is very high. The smooth optimized trajectory also shows the validity of the optimization method with high-order model. Compared with other methods with the 3-DOF model, such as the chaotic particle swarm optimization (CPSO) method and the Gauss pseudospectral method, the trajectory optimization with the high-order model of the parafoil delivery system demonstrates some advan-

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.3 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

3

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80 81

15

Fig. 1. Parafoil delivery system.

16

82 83

17 18 19

84

Table 1 Parameters of the parafoil delivery system.

85

20

Parameter

Value (unit)

86

21

Span Chord Aspect ratio Area of canopy Rope length Mass of canopy Mass of payload Characteristic area of drag of payload Average horizontal velocity Average vertical velocity

10.5 (m) 3.1 (m) 3.4 33 (m2 ) 6.8 (m) 2.1 (kg) 80 (kg) 0.6 (m2 ) 13.8 (m/s) 4.6 (m/s)

87

22 23 24 25 26 27 28 29

88 89 90 91 92 93 94 95 96

30 31 32 33

97

tages. By considering the complicated dynamic constraints, unlike the idealized conditions used with the 3-DOF model, the optimized trajectory is much more practical.

98 99 100

34 35

101

2. Preliminaries

102

36 37 38 39 40 41 42 43 44 45

A parafoil delivery system is a special UAV. The payload and the parafoil canopy are connected by suspension lines [35]. We design a parafoil delivery system whose structure is shown in Fig. 1. In this paper, to realize accurate trajectory optimization, it is necessary to establish a model. Both the 3-DOF model and 6-DOF dynamic model are presented in this section. The parameters of the parafoil delivery system are set as the actual system that is applied in our experiments. The parameters are shown in Table 1.

46 47

2.1. Traditional 3-DOF model

48 49 50 51 52 53 54 55 56 57 58 59 60 61

In the 3-DOF model, the flight velocity and flight direction are primarily considered. It can be expressed as:

⎧ x˙ = v x cos(ψ) + v wx ⎪ ⎪ ⎨˙ y = v y cos(ψ) + v w y ˙ z = vz ⎪ ⎪ ⎩ ˙ ψ =u

(1)

where [x, y ] denotes the horizontal position of the system; z denotes the flight altitude; ψ denotes the yaw angle; u denotes the control quantity of the 3-DOF model; v x , v y denotes the flight velocity in x and y coordinates, respectively; v z denotes the vertical velocity; v wx and v w y denotes the wind velocity in x and y coordinates, respectively.

62 63

2.2. Complicated dynamic constraints of parafoil delivery system

64 65 66

In the 3-DOF model, as shown in Eq. (1), the control quantity directly denotes the angular velocity of the yaw angle. Further-

103 104 105 106 107 108

Fig. 2. Flight trajectory with a high windy environment.

109 110 111

more, wind disturbance is regarded as a component of the flight velocity. In addition, the flight velocity is constant in the 3-DOF model, and the relationship between the control quantity and the angular velocity of the yaw angle is linear. However, in the actual system, both the horizontal and vertical flight velocities will be changed by the horizontal control quantity (flap deflection), the external wind disturbance and other environmental uncertainties. The relationship between the flap deflection and the angular velocity of the yaw angle is also nonlinear. As shown in Fig. 2–4, the complicated dynamic constraints can be observed in the results of the actual flight experiment. Fig. 2 and Fig. 3 present the flight trajectories with a constant horizontal control quantity on the left side. In these two experiments, the wind disturbances are different. Fig. 4 presents the flight trajectory with the trajectory tracking of a straight line. From Fig. 2 and Fig. 3, the effect of wind disturbance on the trajectory and the flight velocity can be clearly observed. A strong wind disturbance will have a large influence on the flight trajectory, as shown in Fig. 2(b). Meanwhile, as shown in Fig. 2(c)-(d) and Fig. 3(c)-(d), both the horizontal and vertical velocities will change with wind disturbance.

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

4

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.4 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15

81

Fig. 5. Coordinate frames and the segmenting of parafoil.

82

16 17 18 19 20 21 22 23 24 25 26

Fig. 3. Flight trajectory with a low windy environment.

27 28 29 30 31 32 33

⎡

34 36 37 38 40 41 42

45 46 47 48 49 50 51

Fig. 4. Flight trajectory with trajectory tracking.

53

57 58 59 60 61

Furthermore, as shown in Fig. 4(b), the parafoil delivery system has large inertia. With such a large inertia, it requires nearly 40 s to stabilize the flight trajectory and track the reference path. All of these dynamic characteristics of the system form the dynamic constraints in the flight trajectory. Meanwhile, the other constraints, including the initial (terminal) constraint, control constraint, and terrain avoidance, are all essential for trajectory optimization of the parafoil delivery system.

62 63

3. Dynamic 6-DOF model of parafoil delivery system

64 65 66

Recently, various high-order models for the parafoil delivery system have been proposed, such as a 6-DOF model from Mor-

83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100

(2)

101 102 103

cθ cψ T e −s = ⎣ s φ s θ c ψ − c φ s ψ cφ s θ c ψ + s φ s ψ

44

56

⎤

⎡

43

55

⎡

where V x , V y , V z denote the flight velocity in the x, y, z axes of the parafoil coordinate system, respectively; T e−s denotes the transformation matrix from geodetic coordinate to parafoil coordinate:

39

54

⎤

x˙ Vx ⎣ y˙ ⎦ = T eT−s ⎣ V y ⎦ z˙ Vz

35

52

taloni [36], an 8-DOF model from Zhu [37] and the 9-DOF models from Prakash [38], Slegers [39] and Mooij [40]. In this paper, on one hand, the complicated dynamic constraints must be considered, therefore the high-order model needs to be applied. On the other hand, the realizability and convergence of the algorithm are also important. The trajectory must converge to the optimal solution in finite time. Synthesizing the above conditions, a complete 6-DOF model is applied for trajectory optimization. In the model, the payload and the parafoil canopy are regarded as an entirety, so that the aerodynamic characteristics of the parafoil can be considered. The coordinate frames are shown in Fig. 5(a). The parafoil coordinate is O s xs y s z s and the geodetic coordinate is O e xe y e ze . O s is also the mass center of the system. First of all, the flight velocity of the system in geodetic coordinate can be expressed by the velocity in parafoil coordinate:

c θ sψ sφ sθ sψ + c φ c ψ c φ sθ sψ − sφ c ψ

−sθ

⎤

sφ c θ ⎦ cφ cθ

104 105 106 107 108 109

(3)

110 111 112

where φ , θ and ψ denote the roll, pitch and yaw angle in geodetic coordinate, respectively. Meanwhile, define sin(x) ≡ s x and cos(x) ≡ c x . The Euler angular velocity in geodetic coordinate can be expressed as:

113

⎡

118

⎤

⎤⎡

⎡

⎤

1 sφ t θ cφ tθ φ˙ Wx ⎣ θ˙ ⎦ = ⎣ 0 cφ −sφ ⎦ ⎣ W y ⎦ Wz 0 sφ c θ−1 c φ c θ−1 ψ˙

114 115 116 117 119

(4)

120 121 122

where W x , W y and W z denote the angular velocity in parafoil coordinate. Define tan(x) ≡ t x . Then, the 6-DOF dynamic model can be outlined as:

123

⎡

126

V˙ x ⎢ V˙ y ⎢ ⎢ V˙ z ⎢ ⎢W ˙ ⎢ x ⎣W ˙y ˙z W

⎤ ⎥ ⎥ ⎥ ⎥ = A 11 ⎥ A 21 ⎥ ⎦

124 125 127

A 12 A 22

− 1

F M

128

(5)

129 130 131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.5 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5

where F denotes all the forces; M denotes all the moments of the forces; A 11 denotes the apparent mass of the parafoil and the real mass of the system; A 22 denotes the rotational inertia; and T A 12 = − A 21 denotes the coupling term. In Eq. (5), the mass matrix can be expressed as:

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

A 11 = ms I 3×3 + ma A 22 = I s + I a −

(6)

× L× s− p ma L s− p

(7)

A 12 = −ma L × s− p

⎡

⎤

− z s− p

0

⎣ z s− p L× s− p = − y s− p

(8)

y s− p −xs− p ⎦ 0

0 xs− p

(9)

where ms denotes the real mass of the total system; ma denotes the apparent mass [41]; I 3 X3 denotes the unit matrix; I s and I a denote the rotational inertia of the real mass and apparent mass, respectively; L × s− p denotes the third order rotation matrix from the mass center O s to the apparent mass center of the parafoil canopy  T y s− p z s− p P , l× denotes the vector from the mass s− p = xs− p center O s to the apparent mass center of the parafoil canopy P . Meanwhile, the force and the moment of the force can be expressed as follows:

F = Fg +

s F aero

+

w F aero

+

r F extr

+

a F extr

(10)

s w r M = M aero + M aero + M extr + M aextr

(11)

28 29 30 31 32 33 34 35 36 37

where F g denotes the gravity of the system; subscript aero denotes the parameters of the aerodynamic force; superscripts s and w denote the parameters of the parafoil and payload, respectively; r and a denote the coupling factors of the real mass and apparent mass, respectively; the subscript extr denotes the coupling force of the system. Furthermore, as shown in Fig. 5(b), in the calculation of the aerodynamic force of the parafoil, the canopy is divided into eight parts. The total aerodynamic force of the parafoil is the sum of these eight parts:

38 39 40

s F aero =

8 

42 44 45 46 47 48 49 50 51

s M aero =

8  i =1

L× T ( F + F di ) s−i i −s li

54

(13)

where F l , F d denote the lift force and drag force, respectively; i denotes the serial number of each segment of canopy; and T i −s denotes the transformation matrix from each segment to O s ; L × i −s

is defined as same as L × s− p . In Eqs. (10)–(13), the aerodynamic force and moment can be expressed as:



w F aero = −C d 0.5ρ S w

57 58 59

w M aero

=



 2

V x2 + V y2 + V z

w L× s− w F aero



F li = ki C li 0.5ρ S i



 2

V x2 + V z

60 61 62 63 64 65 66

C li = f li (α , ul , u r )

(18)

C di = f di (α , ul , u r )

(19)

F di = −C di 0.5ρ S i



⎡

+

V y2

⎡

Vz ⎣ 0 ⎦ −V x

+

V z2 ⎣

⎤

Vx ⎣ Vy ⎦ Vz

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

where α denotes the angle of the attack; ul and u r denote the flap deflection on the left and right side, respectively. ul and u r are also the control quantity of the system. In this way, combining the above formulas, the 6-DOF dynamic model of the parafoil delivery system is established. The complicated dynamic constraints, especially the apparent mass and the aerodynamic characteristics of the canopy, are all considered.

88 89 90 91 92 93 94 95 96

4. Trajectory optimization based on Gauss pseudospectral method

97 98 99

In recent years, considering that the Gauss pseudospectral method can achieve a global optimal control objective, it has been widely used in trajectory planning for many kinds of vehicles [42, 43], such as differentially flat hypersonic vehicle systems [44] and a ramjet-powered vehicle [45]. In this paper, trajectory optimization is realized by combining the 6-DOF dynamic model and the Gauss pseudospectral method. In this section, first, the design of the optimization method and the objective function are presented. Then, the principle of the Gauss pseudospectral method is explained.

100 101 102 103 104 105 106 107 108 109

(16)

⎤

Vx Vy ⎦ Vz

111 112

To achieve trajectory optimization with the high-order model, the toughest and most complicated work is the design of the optimization method, especially the selection of the system states and their updating equations. The states must satisfy the requirements of the flight mission, and the updating method must fit the optimization algorithm.

113 114 115 116 117 118 119

(14) (15)

⎤

⎡

V x2

67

110

55 56

area of the drag force of the payload; i denotes the serial number of each segment; ki denotes the product factor of each segment; S i denotes the characteristic area of the segment; C l and C d denote the parameters of the lift and drag force, respectively. Furthermore, in the model of parafoil delivery system, flight direction control is achieved by the deflection of the trailing edge (flap) of the parafoil canopy. It is realized by pulling down the control rope attached to the rudder of the system. The changing shape of the canopy will lead to complex changing of the aerodynamic force of the parafoil. Thus, its flight direction can be controlled. The deflection of the trailing edge, so called flap deflection, is also the only control quantity of the system. In this paper, C l and C d can be deduced from the attack angle, the control quantity (deflection angle of the flap), and other influences [1]

4.1. Design of the trajectory optimization

52 53

(12)

i =1

41 43

T i −s ( F li + F di )

5

4.1.1. States of the system Utilizing the 6-DOF model, to consider the complicated dynamic constraints, the trajectory optimization method introduces twelve states of the system. These states can be expressed as:



X= x

y

z

φ θ

ψ

Vx

Vy

Vz

Wx

Wy

Wz

T

where C d denotes the drag coefficient of the payload; ρ denotes the real-time air density; S w denotes the characteristic

121 122 123 124 125

(20) (17)

120

In Eq. (20), to achieve accurate homing control, as a basis, the flight position [x, y , z] in the geodetic coordinate system must be involved in the optimized states. It presents the real-time location of the system. Similarly, as the basic states of the system, the roll, pitch and yaw angles [φ, θ, ψ] in the geodetic coordinate

126 127 128 129 130 131 132

JID:AESCTE

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.6 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

6

1 2 3 4 5 6 7 8

system are also introduced into the optimization. The yaw angle not only presents the flight direction of the system, but also will be applied to achieve a flared landing. Then, by considering the complicated dynamic constraints of the system, utilizing the 6-DOF dynamic model, the flight velocity [ V x , V y , V z ] and angular velocity [ W x , W y , W z ] in the parafoil coordinate are all involved. These six states are updated by the model equations and applied to calculation of [x, y , z] and [φ, θ, ψ].

9 10 11 12 13 14 15 16

4.1.2. Updating the states In the optimization, all twelve states in Eq. (20) should be considered. Therefore, these states must be updated by their rate of change:



X˙ = x˙

y˙

z˙

φ˙ θ˙ ψ˙

V˙ x

V˙ y

V˙ z

˙x W

˙y W

˙z W

T

17

(21)

18 19 20

In this section, updating the states can be expressed as:

21 22

X (t + 1) = X (t ) + X˙ (t + 1) T

(22)

23 24 25 26 27 28 29 30 31 32 33 34 35

where T denotes the sampling time. First, the states in geodetic coordinate can be deduced from the states in parafoil coordinate. From Eq. (2), [˙x, y˙ , z˙ ] can be updated by

⎧ x˙ (t + 1) = c θ (t ) c ψ(t ) V x (t ) + (sφ(t ) sθ (t ) c ψ(t ) − c φ(t ) sψ(t ) ) V y (t ) ⎪ ⎪ ⎪ ⎪ cφ(t ) sθ (t ) c ψ(t ) + sφ(t ) sψ(t ) ) V z (t ) +( ⎪ ⎪ ⎨ y˙ (t + 1) = c θ (t ) sψ(t ) V x (t ) + (sφ(t ) sθ (t ) sψ(t ) + c φ(t ) c ψ(t ) ) V y (t ) ⎪ ⎪ +( ⎪ c φ(t ) sθ (t ) sψ(t ) − sφ(t ) c ψ(t ) ) V z (t ) ⎪ ⎪ ⎪ ⎩˙ z(t + 1) = −sθ (t ) V x (t ) + sφ(t ) c θ (t ) V y (t ) + c φ(t ) c θ (t ) V z (t )

36

(23)

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

In Eq. (23), it can be seen that updating of [˙x, y˙ , z˙ ] requires [ V x , V y , V z ] and [φ, θ, ψ]. ˙ can be deduced from Eq. (4): ˙ θ˙ , ψ] Then, [φ,

⎧ ˙ t + 1) = W x (t ) + sφ(t ) t θ (t ) W y (t ) + c φ(t ) t θ (t ) W z (t ) ⎨ φ( ˙θ(t + 1) = c φ(t ) W y (t ) − sφ(t ) W y (t ) ⎩ ˙ ψ(t + 1) = sφ(t ) c θ−(1t ) W y (t ) + c φ(t ) c θ−(1t ) W z (t )

(24)

In Eq. (24), the updating of [ζ, θ, ψ] requires [ W x , W y , W z ] and [φ, θ, ψ]. In comparison with the simple equations in Eq. (1), it can be seen that, in Eq. (23) and Eq. (24), the calculation method associated with these twelve states is very complex. However, even in this condition, Eq. (23) and Eq. (24) are still relatively simple, compared with the updating equations of [ V x , V y , V z ] and [ W x , W y , W z ]. Combining Eq. (5)-Eq. (17), it presents the most important equations of the 6-DOF model. For example, it can be seen that, the calculation method of the aerodynamic force of the parafoil canopy consists of eight parts, F li and F di . Every part is a matrix of three rows and one column, which is calculated by [ V x , V y , V z ]. Meanwhile, the calculation of C l and C d includes more states of the system, as shown in Eq. (18)–(19). Therefore, to update states X , all time-varying variables in X˙ must be expressed by the components in Eq. (21). The other constant parameters will be expressed as a constant. In this way, by analyzing the sys˙ x, W ˙ y, W ˙ z can be expressed tem model in detail, V˙ x , V˙ y , V˙ z , W as:

⎧ V˙ x (t + 1) = f V x (t ) = 0.0001475W z (t )(100.4W x (t ) ⎪ ⎪ ⎪ ⎪ − ⎪ ⎪ 64.26W z (t ) − 9.623sθ (t ) + 0.9821V y (t ) W z (t ) ⎪ ⎪ −0.9821V z (t ) Wy (t ) ⎪  ⎪ ⎪ ⎪ + 0.09705W z (t ) V y (t ) + 0.1345W x (t ) · · · · · · ⎪ ⎪ ⎪ ⎨ V˙ (t + 1) = f (t ) = 8.919c s y Vy θ (t ) φ(t ) − 0.9101V x (t ) W z (t ) + 0 . 9101V ( t ) W ( t ) − ( 2 . 628e − 5) · V x (t ) W x (t ) ⎪ z x ⎪ ⎪ ⎪ ⎪ +0.0001187W x (t ) V x (t ) − 0.0001187V y (t ) W y (t ) · · · · · · ⎪ ⎪ ⎪ ⎪ V˙ z (t + 1) = f V z (t ) = 0.7177V x (t ) W y (t ) − 0.7177V y (t ) W x (t ) ⎪ ⎪   ⎪ ⎪ −0.07093W x (t ) V y (t ) + 0.1345W x (t ) + 0.1345W x (t ) ⎪ ⎪   ⎩ +0.01376W y (t ) V x (t ) − 3.13W y (t ) · · · · · · (25)

⎛

67 68 69 70 71 72 73 74 75 76 77 78 79 80

˙ x (t + 1) = f W x (t ) = 0.003314V x (t ) W z (t ) W ⎜ −0.007341cθ (t ) sφ(t ) − 0.003314V z (t ) W x (t ) ⎜ ⎜ +0.002173V x (t ) W x (t ) − 0.009813W x (t ) V x (t ) ⎜ ⎜ +0.009813V y (t ) W y (t ) · · · · · · ⎜

81 82 83 84

⎜ ˙ ⎜ W y (t + 1) = f W y (t )  ⎜ ⎜ = 0.002505W z (t ) 100.4W x (t ) − 64.26W z (t ) ⎜ ⎜ −0.1216sθ (t ) + 0.01475V y (t ) W z (t ) − 0.01475V z (t ) W y (t ) ⎜ ⎜ +0.001458W z (t ) V y (t ) + 0.1345W x (t ) · · · · · · ⎜ ⎜W ⎜ ˙ z (t + 1) = f W z (t ) = 0.00438c θ (t ) sφ(t ) ⎜ −0.001977V x (t ) W z (t ) + 0.001977V z (t ) W x (t ) ⎜ ⎝ −0.01363V x (t ) W x (t ) +0.06158W x (t ) V x (t ) − 0.06158V y (t ) W y (t ) · · · · · · (26)

85 86 87 88 89 90 91 92 93 94 95 96

In this way, unlike the 3-DOF model, the states of the system, especially the flight velocity and the angular velocity, are deduced from the dynamic model with Eq. (21)–(26). Different from the optimized trajectory with the simple 3-DOF model, the trajectory in this method will be formed by considering the dynamic characteristics of the system model. Therefore, the dynamic constraints can be considered in the trajectory optimization. This is one of the highlights of this paper. Note that, the entire expressions of the flight velocity and angular velocity are much bigger than Eq. (25)–(26). They are composed of all the states in Eq. (21) and the control quantity. Each component of Eq. (25) and Eq. (26) includes at least 328 state variable terms. This complicated expression of the states dramatically challenges the computing power and convergence of the optimization method.

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112

4.1.3. Multiple constraints In this paper, we not only consider the dynamic constraints of the parafoil delivery system, but also involve the environmental constraints and multiple control objectives, such as the terrain avoidance, optimal control quantity, and wind disturbance. 1. Initial constraint and terminal constraint Define t 0 as the initial time and t f as the final time. Then, the initial constraint of the states is set first. It can be expressed as:

⎧ x(t 0 ) = x0 ⎪ ⎪ ⎪ y (t 0 ) = y 0 ⎪ ⎪ ⎨ z(t 0 ) = z0 , ζ (t 0 ) = φ0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ θ(t 0 ) = θ0 ψ(t 0 ) = ψ0

⎧ V x (t 0 ) = V x0 ⎪ ⎪ ⎪ V (t ) = V y 0 ⎪ ⎪ ⎨ y 0 V z (t 0 ) = V z 0 W ⎪ x (t 0 ) = W x0 ⎪ ⎪ ⎪ W ⎪ ⎩ y (t 0 ) = W y 0 W z (t 0 ) = W z 0

113 114 115 116 117 118 119 120 121 122 123 124 125

(27)

126 127 128 129 130

The initial constraint actually denotes the initial states of the system. All the initial constraints should be set before optimization.

131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.7 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

Meanwhile, the terminal constraint is designed as:

1 2 3 4 5 6

⎧ ⎪ ⎪ x(t f ) = 0 ⎨ y (t f ) = 0 ⎪ z(t f ) = 0 ⎪ ⎩ ψ(t f ) = −ψwdir

(28)

7

Furthermore, considering that the range of the control quantity is [−10, 10], this design also can reduce the oscillation of the control quantity and the flight path. In this way, the optimization design is achieved. Then, in the next section, the principle of the Gauss pseudospectral method is presented.

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

where ψwdir denotes the direction of the wind disturbance. In this part, only the specific states need to be defined in the terminal constraint. First, the landing position is designated as the origin of the geodetic coordinate, therefore, the terminal value of the position is [x, y , z] = [0, 0, 0]. By setting the terminal value of the landing position, the objective of landing precision can be ensured. Moreover, in the terminal constraint, the landing direction should be opposite to the direction of the wind disturbance. In this way, the horizontal and vertical flight velocities will be reduced to a great extent, which can help protect the payload. Therefore, as shown in Eq. (28), the yaw angle of the system should be opposite to the wind direction ψwdir . 2. Control constraint In this part, unlike the turning angle of the 3-DOF model, the control quantity (the deflection angle of the flap) is the actual control quantity. It has two sides, that is, left and right. Define u max = 10. Then,



0 ≤ ul ≤ u max 0 ≤ u r ≤ u max

(29)

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

3. Real-time path constraint In this paper, considering that operation of the parafoil delivery system includes mountainous areas and other complex terrain environments, the path constraint is also introduced into the system. The mountain disturbance denotes an unreachable position, and the mountain is designated as a circle area. Then,

  (x − x p , y − y p ) < R p 2     π z = h p sin R p − ( x − x p )2 + ( y − y p )2 δ

55 56

(30) (31)

where p denotes the serial number of the mountain; [x p , y p ] denotes the center of a mountain; h p denotes the altitude of the mountain; δ denotes the smooth coefficient of the mountain. The path constraint is expressed as

     x(t ) − x p , y (t ) − y p  ≥ R p − δ arcsin z(t ) 2

π

53 54

4.2. The principle of the Gauss pseudospectral method

59 60 61 62

hp

(32)

To prevent collision with the mountain, with the path constraint, the control quantity cannot allow the system to break into the mountain area. 4. Objective function The main objectives, including the landing precision, landing direction, and terrain avoidance, have already been considered in the above constraints. Therefore, in this section, the energy management is considered

63

66

71 72

75

u 2 dt

J= t0

76 77 78

t f





J = Φ x(t f ), t f +



 g x(τ ), u (τ ) dτ

79

(34)

82

It must be subject to the dynamic equations

dx dt



= f x(t ), u (t ), t

83 84



(35)

85 86 87

The boundary equality constraints must satisfy

  φmin ≤ φ x(t 0 ), t 0 , x(t f ) ≤ φmax

88

(36)



89 90

The inequality path constraints must satisfy



80 81

t0

91 92

C x(τ ), u (τ ), t 0 , t f ≤ 0

(37)

4.2.1. Discretization of the state variable and control variable First, the discrete point of the Gauss pseudospectral method is distributed in [−1, 1]. Therefore, the real time [t 0 , t f ] needs to be transformed to [−1, 1]. By involving the time variable τ , we have:

τ=

93

(33)

(38)

t f − t0

d  N

1 2N N

dx K

x2 − 1

N 

97 98 100

x(τ ) ≈ Xˆ (τ ) =

N 

L i (τ ) Xˆ i

105 106 108 110 111 112 113 114 115 116 117 118 119

(40)

121 122 123

τ − τj , i = 0, 1 , . . . , N τ i − τj j =0, j =i

124

(41)

125 126

Meanwhile, the control variable is set at the N points in the range of (−1, 1). It is designated similarly

k =1

104

120

N 

u (τ ) ≈ Uˆ (τ ) =

103

109

i =0

N 

102

107

(39)

In this paper, considering that the proposed method also involves initial and final value constraints, there are two more points. Define τ0 and τ N +1 = τ f as the initial and final time, respectively. We then have [τ0 , τ1 , τ2 , . . . , τk , . . . , τ N , τ N +1 ]. Finally, there are N + 2 points. The state variable will be discretized at the first N + 1 points. It can be expressed as [ Xˆ 0 , . . . , Xˆ N ]. After discretizing, the value at these points denote the actual value. Then, x(τ ) can be obtained by forming Lagrange interpolation polynomials to approximate the state variable

L i (τ ) =

96

101

In the open interval [−1, 1], it takes N points (Legendre-Gauss (LG) points). Define [τ1 , τ2 , . . . , τk , . . . , τ N ]. These points are the zero point of the N order Legendre polynomial. The N order Legendre polynomial can be expressed as:

H N (x) =

95

99

2t − t f − t 0

t f

64 65

70

74

In the Gauss pseudospectral method, the form of the optimal problem can be expressed as follows:

57 58

69

94

where ul and u r denote the control quantity of the left and right side, respectively. Unlike the 3-DOF model, the influence of the control quantity on the yaw angle is not set. The influence will be much more complicated in the 6-DOF model. It will be evaluated and optimized in GPM.

34 35

68

73

7 8

67

127 128 129 130

L k∗ (τ )Uˆ k

(42)

131 132

JID:AESCTE

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.8 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

8

1

67

Table 2 Simulation environment.

2 3

Simulation

4 5

Case 1 Case 2 Case 3

6 7

68 69

Constant wind

Gust disturbance

Terrain (Mountain)

Velocity

Direction

Velocity

Direction

Center

Radius

Altitude

3 m/s 3 m/s 3 m/s

0 (rad) 0 (rad) 0 (rad)

3 m/s 3 m/s 3 m/s

0 (rad) 3.14 (rad) 3.14 (rad)

[500, 1000] [1500, 1500]

350 350

1500 1500

70 71 72 73 74

8 9 10 11

L k∗ (τ ) =

j =1, j =k

12 13

N 

τ − τj , i = 1, 2, . . . , N τk − τ j

(43)

x˙ (τ ) =

where L i and L k∗ denote Lagrange interpolation functions.

14 15 16 17 18 19 20 21 22

4.2.2. Initial and terminal constraints From Eq. (40), it can be observed that the corresponding points of x(τ ) are the first N + 1 points. x(τ ) does not have initial and terminal control quantities. However, in this paper, the optimization problem includes the initial and terminal constraints. The initial constraint is known and set first. Then, the terminal constraint of the system can be expressed as:

23 24 25

Xˆ f = Xˆ 0 +

t f − t0 2

28 29 30

33 34 35

38

Xˆ f = Xˆ 0 +

2

41

(44)



ω¯ k f Xˆ (τ ), Uˆ (τ ), τk , t 0 , t f



(45)

k =1

u (τ0 ) ≈ Uˆ (τ0 ) =

N 

L k∗ (τ )Uˆ k

(46)

k =0

u (τ f ) ≈ Uˆ (τ f ) =

N +1 

L k∗ (τ )Uˆ k

(47)

k =0

42 43 44 45 46 47

4.2.3. Transformation of the equation of motion The equation of motion of the system must be transformed to the derivation operation of the Lagrange interpolation functions. After the derivative of Eq. (40), we have

48 49 50

53 54

L˙ i (τ ) Xˆ i

N 

D ki = L˙ i (τk ) =

63

66

N  j =0, j =i

(τi − τ j )

k = 1 , 2 , . . . , N , i = 0, 1 , . . . , N

(49)

In this way, the derivative of the state variable

64 65

j =0, j =i ,l

l =0

60 62

N 

(τk − τ j )

x˙ (τ ) =

N  i =0

L˙ i (τ ) Xˆ i =

N  i =0

D ki (τk ) Xˆ i

(50)

75 76 77 78 79

In this way, the equation of motion in Eq. (50) can be replaced by Eq. (51).

80

4.2.4. Real-time path constraint and discretization of the objective function By the discretization of the real-time path constraint at point τk , it can be obtained that

83

C p ( Xˆ k , Uˆ k , τk , t 0 , t f ) ≤ 0,

(52)

88

Then, the objective function also needs to be dispersed at point

90

p = 1, 2, . . .

τk

81 82 84 85 86 87 89 91 92

+

t f − t0 2

N 

93 94

ω¯ k L ( Xˆ k , Uˆ k , τk , t 0 , t f )

(53)

k =1

With the above transformation, the problem of trajectory optimization with complicated dynamic constraints and external environmental influence can be regarded as a series of nonlinear parameter optimization problems. Meanwhile, the sequential quadratic programming algorithm is applied for the optimization. The optimal trajectory can be obtained by the discrete calculation. In the process of solving all these constraints with the Gauss pseudospectral method, a small number of points are first calculated to build the draft of the trajectory. Then, the other points are inserted into these calculated points to smooth the trajectory by the interpolation method. In this way, the precision of the optimized trajectory can be guaranteed while maintaining a sufficiently fast computation speed. 5. Simulation experiment 5.1. Simulation environment

A differential matrix can be obtained by the derivative of the Lagrange interpolation polynomials at point τk

59 61

(51)

1. Simulation with different initial positions: In this part, three different initial positions are applied. The initial positions are [1000, 3000, 1500], [2500, 2500, 1500] and [−1500, 1500, 1200]. The results will illustrate the convergence of the optimization method with different initial values. The landing position is [0, 0, 0]. There are two mountains on the map. For both mountains, the altitude is 3000 m and the radius is 400 m. The centers are [500, 1000] and [1500, 1500]. In the simulation, there are two parts of the wind disturbance: the constant wind and gust. The known constant wind of 3 m/s from the negative direction of the X coordinate relative to the geodetic coordinate system was applied throughout the entire flight process. Meanwhile, to test the optimization ability with environmental uncertainty, a gusty wind model was also introduced into the simulation. It is unknown and will directly change the flight velocity of the system. It is applied

i =0

56 58

f ( Xˆ k , Uˆ k , τk , t 0 , t f )

In this section, the simulation environment is presented as Table 2. The simulation is carried out with three parts:

55 57

2

(48)

x˙ (τ ) =

51 52

N 

t f − t0

J ∗ = Φ( Xˆ 0 , Xˆ f , t 0 , t f )

¯ k denotes the integral weight factor. where ω Then, the initial and terminal control quantities can be expressed as:

39 40

f x(τ ), u (τ ), τ dτ

N t f − t0 

36 37



After discretizing, we have

31 32



−1

26 27

1

Meanwhile, after discretization of the equation of motion, it can be obtained that

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.9 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17 18 19 20

83

Fig. 7. Control quantity with different initial positions.

Fig. 6. Horizontal flight trajectory with different initial positions. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

84 85 86 87

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

to verify the capability of the wind disturbance resistance in trajectory optimization based on the Gauss pseudospectral method. In Case 1, a gust of 3 m/s from the negative direction of the X coordinate relative to the geodetic coordinate system was applied at time 150 s from the exit, for a total time of 15 s. In Case 2, the gust of 3 m/s from the positive direction of the X coordinate was applied at time 60 s from the exit, for a total time of 15 s. In Case 3, the gust of 3 m/s from the negative direction of the X coordinate relative to the geodetic coordinate system was applied at time 100 s from the exit, for a total time of 15 s.

88 89 90 91 92 93 94 95 96 97

2. Simulation with different terrain environments: The simulation was carried out with different terrain disturbances. The initial position is [1000, 3000, 1500]. The wind disturbance is the same as Case 2. This part is applied to verify the capacity of the terrain avoidance of the proposed method. It will also present the influence of the terrain disturbance on the dynamic constraints.

98

3. Comparison experiment: In this section, the proposed method is compared with two other optimization algorithms. These two optimization algorithms, CPSO and the Gauss pseudospectral method, are designed with the 3-DOF model introduced in Section 2.1. The comparison with the CPSO presents the advantage of the calculation ability of GPM. Meanwhile, the comparison with the GPM using the 3-DOF model presents the advantage and improvement of the proposed method by considering the complicated dynamic constraints of the parafoil delivery system.

106

99 100 101 102 103 104 105 107 108 109 110 111 112 113 114 115 116

50 51

5.2. Simulation with different initial positions

117

Fig. 8. 3D flight trajectory with different initial positions.

118

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

First, in the simulation experiment, the airdrop experiment with different initial positions was carried out. This section focuses on the influence of the initial value in the trajectory optimization. The dynamic disturbances, including the wind disturbance and terrain disturbance, are all considered. The results are presented in Fig. 6–10. Fig. 6–9 presents the results of the Case 1 and Case 2. Fig. 6 shows the horizontal flight trajectory. Fig. 7 shows the control quantity. Fig. 8 shows the 3D flight trajectory. Fig. 9 shows the angular velocities W x , W y and W z . Fig. 10 presents the horizontal flight trajectory and control quantity of Case 3. From Fig. 6 and Fig. 8, it can be observed that the trajectory is optimized successfully with the proposed method. Meanwhile, all the dynamic constraints are satisfied under wind and gust disturbances. First, considering that the state of the final position is

regarded as the terminal constraint in the Gauss pseudospectral method, the landing error is 0 m in both cases. Fig. 6(a) shows that the terrain avoidance is completed as well. In the magnified part of Fig. 6, it can be observed that the landing direction is opposite to the constant wind direction to reduce the flight velocity in the landing stage. Most importantly, the optimized trajectory is smooth with the complicated dynamic constraints of the parafoil delivery system. Therefore, it is practicable and easy to track in the trajectory tracking system. In Fig. 7, u presents the control quantity - deflection angle of the flap. It is expressed as:



ul = u , u r = 0 u l = 0, u r = − u

if u > 0 if u ≤ 0

119 120 121 122 123 124 125 126 127 128 129 130 131

(54)

132

JID:AESCTE

10

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.10 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101 102

36 37

Fig. 9. Angular velocity with different initial positions.

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

103 104

38

Considering that the constant wind disturbance is from the negative direction of the X coordinate, the flap deflection for both cases is negative without the gust disturbance. This allows the system to turn right and resist the constant wind disturbance. Moreover, at 150 s in Case 1 and 60 s in Case 2, it can be observed that the control quantity has a quick reaction to the unknown gust disturbance, which proves that the proposed method has the ability to resist external variable disturbances. Furthermore, by considering the complicated dynamic constraints, unlike optimization with the 3-DOF model, the flight time is not constant. It is variable based on different control conditions. From the results, the flight time is 298 s for Case 1 and 304 s for Case 2. The different flight time is due to the increased drag associated with the flap deflection. Initially, for deflections less than approximately 15%, the deflection will create lift with a small rolling moment, and the payload suspended below the parafoil will create a counter roll. However, with deflections greater than approximately 15%, the deflection will dominate the aerodynamics of the wing section and result in a yaw, leading to a combined roll-yaw in the direction of the deflected wing section. With the consideration of the constraints of the angular velocity in 6-DOF model, the optimized trajectory with the proposed method is more practical. Fig. 9 presents the angular velocity of the system in parafoil coordinate. It can be observed that, considering the complicated dynamic constraints of the system, the attitude information and the angular velocity show fluctuations under the coupling influence of the constant wind and gust. In Case 1, combining the constant wind and the gust, the deviations of the angular veloc-

ity are produced by an external wind disturbance whose velocity is 6 m/s from the negative direction of X coordinate. The smaller deviations at 120 s in Case 2 are caused by the terrain constraints and the oscillation of the control quantity. Meanwhile, Case 3 with another initial position is presented. The release point is [−1500, 1500, 1200], which is upwind of the target. The results are presented as Fig. 10. Fig. 10(a) shows the flight trajectory. Fig. 10(b) is the control quantity. These results show that the control objective of the landing precision is achieved successfully with the upwind release position. The landing direction of the system is opposite to the constant wind direction, and a flared landing can be realized.

105 106 107 108 109 110 111 112 113 114 115 116 117

5.3. Simulation with different terrain environments

118 119

This section focuses on the influence of the external constraints, especially the environmental constraint, of the optimized trajectory. The initial position is set at [1000, 3000, 1500]. The flight trajectory is shown in Fig. 11 and Fig. 13. The control quantity is shown in Fig. 12. In Fig. 11 and Fig. 13, it can be observed that the trajectory can be obtained under both terrain conditions. The influence of the terrain on the dynamic constraints can be observed from the flight trajectory. The control quantity not only can avoid the terrain disturbance, but also can achieve an accurate landing. Meanwhile, unlike the case with the terrain disturbance, without the terrain disturbance, the redundant flight altitude is consumed by changing its flight direction. This also can be seen from the control quan-

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.11 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

11

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82 83

17 18

Fig. 12. Control quantity with different terrain environments.

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100 101

35 36

102

Fig. 10. Simulation case with the upwind release position.

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116 117

51

Fig. 13. 3D flight trajectory with different terrain environments.

52 54 55

Fig. 11. Horizontal flight trajectory with different terrain environments.

5.4. Comparison experiment with CPSO and Gauss pseudospectral method in 3-DOF model

58 59 60 61 62 63 64 65 66

120 121 122

56 57

118 119

53

tity. In the simulation without the terrain disturbance, the control quantity varies at the initial stage to change the flight trajectory. Even in this way, the trajectory is still smooth. Considering that the control quantity is directly obtained by the optimization method, not a stable controller, the proposed method presents a huge improvement. Furthermore, considering that the flap deflection has a complex influence on the aerodynamic force on the parafoil, even with the same initial position, different control quantities lead to different flight times. The flight time is 296 s for the case without terrain disturbance.

In this section, Case 1 is applied for a comparison experiment. There are three methods: GPM (6-DOF), GPM (3-DOF), and CPSO (3-DOF), which are all applied for the experiment. The results are presented in Fig. 14–15. Fig. 14 presents the horizontal flight trajectory. Considering that the meaning and units of the control quantities of the 3-DOF and 6-DOF models are different, Fig. 15 only presents the control quantities of GPM (3-DOF) and CPSO (3-DOF). The control quantity of the optimization method with the 6-DOF model is presented in Fig. 7. The detailed results using these methods are presented in Table 3.

123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.12 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

12

1 2 3 4 5 6 7

67

Table 3 Comparison of experimental results.

68

Method

Flight time (s)

Landing error (meter)

Landing direction (◦ )

69

GPM (6-DOF) GPM (3-DOF) CPSO 1 (3-DOF) CPSO 2 (3-DOF)

298 326 326 326

0 0 54.7 90.3

−180 −180 −167 −152

70

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Fig. 14. Horizontal flight trajectory for comparison experiment.

27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

Fig. 15. Control quantity for comparison experiment.

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

72 73 74

8

26

71

From Fig. 14 and Table 3, it can be observed that the Gauss pseudospectral method has an advantage with respect to the trajectory optimization problem of the parafoil delivery system. It has better landing precision and an accurate landing direction. In Table 3, the error of the landing position and direction of GPM is zero with the terminal constraints. For better optimization, the highorder model is applied with GPM. Actually, we also tried the CPSO with the 6-DOF model, but it ran for more than 5 hours and still could not obtain a satisfying trajectory. It also illustrates the advantage of GPM. With the complicated aerodynamic equations in the 6-DOF model, it is clear that the computation time will be longer than the optimization case with the 3-DOF model. However, by employing a professional computing device, it is reasonable to achieve the proposed method for the in-flight trajectory optimization. Meanwhile, from Fig. 14, the advantage of the 6-DOF model is presented. First, with the constant vertical flight velocity in the 3-DOF model, the flight time of the 3-DOF trajectory is also constant, as shown in Table 3. However, in the flight process, uncertain disturbances, such as the gust, will influence the flight trajectory directly and produce errors compared with the optimized trajec-

tory. Therefore, the flight time is variable, even with the same altitude. In the 6-DOF model, unlike the idealized flight time in the 3-DOF model, this influence can be considered. By considering the complicated dynamic constraints, the optimized trajectory better reflects reality. Second, in the energy management part of GPM with the 3-DOF model, the largest angular velocity of the yaw angle is 0.142 rad/s, and lasts for approximately 20 s. In this situation, it is difficult to separate this part into the straight line for the trajectory tracking. However, with the 6-DOF model, the control quantity is the actual deflection angle of the flap. The angular velocity of the yaw angle is calculated by the aerodynamic equations of the system. The optimized trajectory is thus realistic and realizable. Based on these statements, it can be observed that the proposed method has better landing precision in the airdrop experiment of the parafoil delivery system. The realizability of the proposed method has the advantage of considering the complicated dynamic constraints and shows a huge improvement compared with other methods.

75

6. Conclusion

95

In this paper, a novel trajectory optimization method for parafoil delivery systems was explored. The proposed method was achieved by combining a 6-DOF dynamic model and the Gauss pseudospectral method. The main enhancement of this method is the consideration of complicated dynamic constraints in the trajectory optimization of parafoil delivery system. Moreover, this was the reason for the use of a high-order model in the optimization. By considering the dynamic constraints, the optimized trajectory can be made more realizable and consistent with reality. Furthermore, the results also show that the Gauss pseudospectral method is suitable for this problem. It not only results in faster computation speed, but also can achieve a globally optimal solution. From the results, it can be seen that the design of the optimization method is correct and the optimized trajectory is smooth. Furthermore, it can be observed that the proposed method can satisfy all the constraints and the multiple control objectives under windy environment, including the dynamic constraints of the system, terrain avoidance and landing precision. Compared with other methods, it presents a huge improvement for the specific purpose of a parafoil delivery system that requires high landing precision and good reliability.

97

Declaration of competing interest The authors declared that they have no conflicts of interest to this work. References [1] O.A. Yakimenko, Precision Aerial Delivery Systems: Modeling, Dynamics, and Control, AIAA Press, 2015. [2] J. Stein, C. Madsen, A. Strahan, An overview of the guided parafoil system derived from X-38 experience, in: 18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2005, p. 1652. [3] T. Jann, Advanced features for autonomous parafoil guidance, navigation and control, in: 18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2005, p. 1642.

76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 96 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105631 /FLA

[m5G; v1.261; Prn:13/12/2019; 12:34] P.13 (1-13)

H. Sun et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

[4] M.R. Cacan, M. Costello, M. Ward, et al., Human-in-the-loop control of guided airdrop systems, Aerosp. Sci. Technol. 84 (2019) 1141–1149. [5] M. Ghoreyshi, K. Bergeron, A. Jirasek, et al., Computational aerodynamic modeling for flight dynamics simulation of ram-air parachutes, Aerosp. Sci. Technol. 54 (2016) 286–301. [6] X. Xue, Y. Nakamura, K. Mori, et al., Numerical investigation of effects of angleof-attack on a parachute-like two-body system, Aerosp. Sci. Technol. 69 (2017) 370–386. [7] B.S. Chiel, C. Dever, Autonomous parafoil guidance in high winds, J. Guid. Control Dyn. 38 (5) (2014) 963–969. [8] B. Le Floch, J. How, L. Breger, et al., Trajectory planning for autonomous parafoils in complex terrain, in: 24th AIAA Aerodynamic Decelerator Systems Technology Conference, 2017, p. 3220. [9] D. Carter, S. George, P. Hattis, et al., Autonomous large parafoil guidance, navigation, and control system design status, in: 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2007, 2007, p. 2514. [10] A.J. Calise, D. Preston, Swarming/flocking and collision avoidance for mass airdrop of autonomous guided parafoils, J. Guid. Control Dyn. 31 (4) (2008) 1123–1132. [11] T. Jann, Aerodynamic model identification and GNC design for the parafoil-load system ALEX, in: 16th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2001, 2015. [12] H. Sun, Q. Sun, S. Luo, et al., In-flight compound homing methodology of parafoil delivery systems under multiple constraints, Aerosp. Sci. Technol. 79 (2018) 85–104. [13] A. Rosich, P. Gurfil, Coupling in-flight trajectory planning and flocking for multiple autonomous parafoils, Proc. Inst. Mech. Eng., Part G, J. Aerosp. Eng. 226 (6) (2012) 691–720. [14] L. Fowler, J. Rogers, Bezier curve path planning for parafoil terminal guidance, J. Aerosp. Inform. Syst. 11 (5) (2014) 300–315. [15] M. Ward, A. Gavrilovski, M. Costello, Performance of an autonomous guided airdrop system with glide slope control, in: 21st AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2011, p. 2572. [16] O.A. Yakimenko, I.I. Kaminer, Optimal control and practical GNC algorithms for autonomous low- and high-glide-ratio delivery systems, IFAC Proc. Vol. 37 (8) (2004) 466–471. [17] C. Madsen, R. Sostaric, C. Cerimele, Flight performance, aerodynamics, and simulation development for the X-38 parafoil test program, in: 17th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2003, p. 2108. [18] A. Strahan, Testing of parafoil autonomous guidance, navigation and control for X-38, in: 17th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2003, p. 2115. [19] B.J. Rademacher, P. Lu, A.L. Strahan, et al., In-flight trajectory planning and guidance for autonomous parafoils, J. Guid. Control Dyn. 32 (6) (2009) 1697–1712. [20] T. Yomchinda, J.F. Horn, J.W. Langelaan, Modified Dubins parameterization for aircraft emergency trajectory planning, Proc. Inst. Mech. Eng., Part G, J. Aerosp. Eng. 231 (2) (2017) 374–393. [21] J.R. Cleminson, Path planning for autonomously-guided parafoils: a dynamic programming approach, in: AIAA Aerodynamic Decelerator Systems (ADS) Conference, 2013, p. 1346. [22] S. Culpepper, M.B. Ward, M. Costello, et al., Adaptive control of damaged parafoils, in: AIAA Aerodynamic Decelerator Systems (ADS) Conference, 2013, 2013, p. 1344. [23] M.R. Cacan, E. Scheuermann, M. Ward, et al., Autonomous airdrop systems employing ground wind measurements for improved landing accuracy, IEEE/ASME Trans. Mechatron. 20 (6) (2015) 3060–3070. [24] J. Tao, Q. Sun, H. Sun, et al., Dynamic modeling and trajectory tracking control of parafoil system in wind environments, IEEE/ASME Trans. Mechatron. 22 (6) (2017) 2736–2745.

13

[25] N. Slegers, M. Costello, Model predictive control of a parafoil and payload system, J. Guid. Control Dyn. 28 (4) (2005) 816–821.

67

[26] B. Luders, A. Ellertson, J.P. How, et al., Wind uncertainty modeling and robust trajectory planning for autonomous parafoils, J. Guid. Control Dyn. 39 (7) (2016) 1614–1630.

69

[27] D. Carter, L. Singh, L. Wholey, et al., Band-limited guidance and control of large parafoils, in: 20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2009, p. 2981. [28] T.Z. Gimadieva, Optimal control of a gliding parachute system, J. Math. Sci. 103 (1) (2001) 54–60. [29] N. Slegers, O. Yakimenko, Optimal control for terminal guidance of autonomous parafoils, in: 20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2009, p. 2958.

68 70 71 72 73 74 75 76 77

[30] A.V.S. Babu, V.C. Suja, C.V. Reddy, Three dimensional trajectory optimization of a homing parafoil, IFAC Proc. Vol. 47 (1) (2014) 847–854.

78

[31] B.D. Luders, I. Sugel, J.P. How, Robust trajectory planning for autonomous parafoils under wind uncertainty, in: AIAA Infotech Aerospace Conference, 2013, p. 4584.

80

[32] J. Rogers, N. Slegers, Robust parafoil terminal guidance using massively parallel processing, J. Guid. Control Dyn. 36 (5) (2013) 1336–1345. [33] N. Slegers, J.D. Rogers, Terminal guidance for complex drop zones using massively parallel processing, in: AIAA Aerodynamic Decelerator Systems (ADS) Conference, 2013, p. 1343.

79 81 82 83 84 85 86

[34] N. Slegers, A. Brown, J. Rogers, Experimental investigation of stochastic parafoil guidance using a graphics processing unit, Control Eng. Pract. 36 (2015) 27–38.

87

[35] J. Wachlin, M. Ward, M. Costello, In-canopy sensors for state estimation of precision guided airdrop systems, Aerosp. Sci. Technol. 90 (2019) 357–367.

89

[36] P. Mortaloni, O. Yakimenko, V. Dobrokhodov, et al., On the development of a six-degree-of-freedom model of a low-aspect-ratio parafoil delivery system, in: 17th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, 2003, p. 2105. [37] E. Zhu, Q. Sun, P. Tan, Z. Chen, X. Kang, Y. He, Modeling of powered parafoil based on Kirchhoff motion equation, Nonlinear Dyn. 79 (1) (2014) 617–629. [38] O. Prakash, N. Ananthkrishnan, Modeling and simulation of 9-DOF parafoilpayload system flight dynamics, in: AIAA Atmospheric Flight Mechanics Conference and Exhibit, 2006, p. 6130.

88 90 91 92 93 94 95 96 97

[39] N. Slegers, M. Costello, Aspects of control for a parafoil and payload system, J. Guid. Control Dyn. 26 (6) (2003) 898–905.

98

[40] E. Mooij, Q. Wijnands, B. Schat, 9 DoF parafoil/payload simulator development and validation, in: AIAA Modeling and Simulation Technologies Conference and Exhibit, 2003, p. 5459.

100

[41] T.M. Barrows, Apparent mass of parafoils with spanwise camber, J. Aircr. 39 (3) (2002) 445–451.

99 101 102 103

[42] G.T. Huntington, A.V. Rao, Optimal reconfiguration of spacecraft formations using the Gauss pseudospectral method, J. Guid. Control Dyn. 31 (3) (2008) 689–698.

104

[43] E. Yong, G. Tang, L. Chen, Rapid trajectory optimization for hypersonic reentry vehicle via Gauss pseudospectral method, J. Astronaut. 29 (6) (2008) 1766–1772.

107

105 106 108

[44] J. Wang, H. Liang, Z. Qi, et al., Mapped Chebyshev pseudospectral methods for optimal trajectory planning of differentially flat hypersonic vehicle systems, Aerosp. Sci. Technol. 89 (2019) 420–430.

109

[45] S. Yang, T. Cui, X. Hao, et al., Trajectory optimization for a ramjet-powered vehicle in ascent phase via the Gauss pseudospectral method, Aerosp. Sci. Technol. 67 (2017) 88–95.

112

110 111 113 114 115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132