Next Article in Journal
An Intelligent Control and a Model Predictive Control for a Single Landing Gear Equipped with a Magnetorheological Damper
Previous Article in Journal
Effects of Injector Configuration on the Detonation Characteristics and Propulsion Performance of Rotating Detonation Engine (RDE)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Coupling Dynamics and Three-Dimensional Trajectory Optimization of an Unmanned Aerial Vehicle Propelled by Electroaerodynamic Thrusters

1
Department of Aerospace Engineering, Harbin Institute of Technology, Harbin 150001, China
2
Aerodynamics Research Institute, Aviation Industry Corporation of China, Ltd. (AVIC), Harbin 150001, China
*
Authors to whom correspondence should be addressed.
Submission received: 28 July 2023 / Revised: 2 November 2023 / Accepted: 5 November 2023 / Published: 10 November 2023

Abstract

:
Electroaerodynamic unmanned aerial vehicles (EAD-UAVs) are innovative UAVs that use high-voltage asymmetric electrodes to ionize air molecules and Coulomb force to push these ions to produce thrust. Unlike fixed-wing and rotor UAVs, EAD-UAVs contain no moving surfaces and have the advantages of very low noise, low mechanical fatigue, and no carbon emissions. This paper proposes an EAD-UAV configuration with an orthogonal arrangement of multiple EAD thrusters to adjust the EAD-UAV attitude and flight trajectory through voltage distribution control alone. Based on a one-dimensional dynamic model of an EAD thruster, the attitude–path coupling dynamics of the EAD-UAV were derived. To achieve EAD-UAV flight control for a specified target, the Bezier shaping approach (BSA) was implemented to realize rapid trajectory optimization considering the coupling dynamic constraints. The numerical simulation results indicate that the BSA can quickly procure an optimized flight trajectory that satisfies the dynamic and boundary constraints. Compared with the Gaussian pseudospectral method (GPM), the BSA changes the optimization index of the objective function by nearly 1.14% but demands only nearly 1.95% of the computational time on average. Hence, the improved integrative Bezier shaping approach (IBSA) can overcome the poor convergence issue of the BSA under the continuous acceleration constraint of multi-target flight trajectories.

1. Introduction

Since the first airplane flight took place in 1903, these vehicles have been controlled by moving surfaces such as flaps, propellers, and turbines. Owing to the advantages of good stealth performance, low maintenance difficulty, and low noise, none-control-surface aircrafts (NCSAs) have received a lot of research [1,2,3]. Unmanned aerial vehicles (UAVs) are increasing rapidly in popularity for use in both military and civil applications [4,5,6,7,8], which makes NCS-UAVs an important subject for research on NCSAs. In 2018, the first unmanned aerial vehicle (UAV) without moving surfaces was successfully tested for the first time [9]; it was propelled by an electroaerodynamic (EAD) thruster. The emergence of EAD-UAVs has provided a huge impetus for the development of NCS-UAVs. EAD-UAVs constitute an innovative UAV concept, using high-voltage asymmetric electrodes to ionize air molecules and Coulomb force to push air ions to produce ionic wind, generating thrust. Unlike fixed-wing and rotor UAVs, there are no moving surfaces, such as propellers or turbines, in EAD-UAVs. As solid-state propulsion UAVs, EAD-UAVs have the advantages of very low noise, low mechanical fatigue, and no carbon emissions [10].
EAD propulsion involves generating propulsive forces in air fluid [11,12], which is an approach to handling and moving fluids without the help of moving surfaces, as illustrated in Figure 1. In 1928, Brown proposed a prototype EAD thruster structure, which uses a stack of asymmetric electrodes to generate directional motion in the air under a high voltage of tens of kilovolts [13]. Because EAD-UAVs can fly ultra-silently by generating ionic wind, this innovative UAV concept has garnered increasing research attention and has been employed in several design proposals recently [14,15]. Scholars at the Central Institute of Technology in Paris, France [16] and Yonsei University in Seoul, Korea [17] have studied complex multi-stage stacked propellers, and their research results have proven the possibility of designing complex multi-stage structures to enhance the performance of plasma propellers. Researchers at Sandia National Laboratory in Albuquerque, New Mexico, tested the effects of AC and DC power supplies and positive and negative polarities on propeller performance [18]. In addition to thrust, the thrust-to-power ratio is an important index for measuring the performance of plasma thrusters. Masuyama et al. [19] studied different factors influencing the thrust-to-power ratio and presented a one-dimensional model of the thrust-to-power ratio of plasma thrusters. Monrolin et al. [14] used another index to analyze the performance of plasma thrusters: thrust density. In addition to changing the parameters of single-stage thrusters to improve their performance, this objective can also be achieved through the series parallel connection of multiple thrusters, similar to the series parallel connection of resistors. Series connection means that multiple thrusters are connected at the end and arranged in the same plane; parallel means that the emitters of multiple thrusters are in a single plane, the receivers are in another plane, and the two planes are parallel. Stuetzer et al. [12] considered the effects of parallel and series connections on ion pumps under one-dimensional conditions. Masuyama et al. [19] studied multistage propellers. Colas et al. [16,20,21] conducted series and classification experiments on plasma thrusters and found that the classification structure could produce a greater ion wind velocity. In addition, the feasibility of using EAD devices to generate sufficient propulsion for UAVs has been discussed and questioned [22]. In recent years, considerable progress has been made in lightweight high-pressure generation technology, and the possibility of using plasma propulsion technology to provide the main thrust in UAVs has gradually increased. In 2016, Wynsberghe and Turak proposed using plasma thrusters in stratospheric floating hot air balloons to generate thrust to drive flight [23]. In 2017, Drew and Pister developed and successfully test flew a micro robot propelled by plasma, with a thrust-to-weight ratio of 10 [24]. In 2018, Xu et al. demonstrated that an EAD propulsion system can maintain powered flight by designing and flying a heavier-than-air airplane propelled by an EAD thruster [9]. The investigated EAD-UAV was a fixed-wing airplane with a 5 m wingspan, which was tested 10 times and obtained steady-level flight. The whole power system, including all of the batteries and a pointedly designed ultralight high-voltage (40 kV) power converter, was attached to the body. The actual flight results indicate that traditionally approved limitations in aspects of the thrust-to-power ratio and thrust density [14,15,22], which were once considered to make EAD thrusters unachievable as an approach to airplane propulsion, have been overcome. This research provides a proof of concept for EAD-UAVs, opening up feasibilities for UAVs and aerodynamic devices that are quieter, mechanically simpler, and have no carbon emissions. The unidirectional flight of EAD-UAVs has been proven to be feasible through actual flight experiments. The flight ability of EAD-UAVs has finally been proved through the efforts of many scholars over the years. In 2021, He et al. used EAD thrusters for both the propulsion and yaw control of a blimp, The feasibility of EAD thrusters for attitude control of lighter-than-air aircraft has also been verified [25]. However, attitude and flight trajectory adjustment of heavier-than-air EAD-UAVs without moving surfaces has not been discussed or investigated in depth thus far.
This paper proposes a new EAD-UAV configuration with an orthogonal arrangement of multiple EAD thrusters that can adjust the attitude and flight trajectory of an EAD-UAV through voltage distribution control alone. As illustrated in Figure 2, six EAD thrusters are arranged orthogonally under the fixed wing, and their thrusts, F1F6, can be controlled independently by adjusting the voltage. Through the combined control of the EAD thrusters, the desired control torque and propulsive force can be generated within a certain range, which is mainly limited by the maximum voltage. The EAD-UAV in this paper is designed by adding several groups of EAD thrusters for attitude control based on the UAV in [9]. In this paper, the influence of structural change on the aerodynamic characteristics of UAVs is not considered. The airframe parameters are the same as those in [9], the aerodynamic parameters only consider the lift coefficient and drag coefficient, and the lateral force coefficient and aerodynamic moment coefficients are set to 0. It is assumed that the EAD-UAV can track the aircraft deflection angle and track the inclination angle during flight, and the fuselage axis direction is always consistent with the flight speed direction, so the lift coefficient and drag coefficient remain unchanged. The lift of the EAD-UAV comes from the installation angle of the wing. The drag coefficient is estimated according to the drag and flight speed given in [9], and the lift coefficient is estimated according to the lift-drag ratio given in [9]. The six EAD thrusters of the EAD-UAV are assumed to have the same thrust at the same voltage, and the structural parameters of the EAD thrusters are consistent with those in [9]. The thrust of thrusters 3, 4, 5, and 6 has two directions, which are positive to the right and upward, and negative for the opposite. The + and − of F3, F4, F5, and F6 represent the direction of thrust. The fuselage parameters of the EAD-UAV and the parameters of the EAD propeller are given in Table A1.
The EAD-UAV designed in this work performs trajectory and attitude coupling control by adjusting the thruster voltage. A trajectory planned without considering the dynamic constraints may not adhere to the voltage limitations of the thruster. Therefore, in this study, a trajectory optimization method considering dynamic constraints was adopted. Trajectory-optimizing algorithms typically use optimal control theory. Popular numerical solutions to optimal control theory include direct methods (“discrete first and then optimization”) and indirect methods (“optimize first and then discrete”) [26], both of which need to be based on dynamics and require considerable amounts of integral calculations. The pseudospectral method is generally used in direct approaches and works well for numerous aerospace problems [27,28,29]. Recently, many improved pseudospectral methods have been proposed. For instance, Rogowski et al. applied the Chebysev pseudospectral method to create the trajectory of a glider in a vertical plane [30]. Further, Pepy et al. introduced an optimal algorithm based on an indirect shooting method [31]. The geometric method was also employed in optimal control theory. Babel et al. applied the shortest path algorithm to deal with network optimization [32]. They divided trajectories in space into multiple line elements, which could act as the UAV trajectory when connected. However, all UAVs’ speeds were assumed to be constant, so the path length was chosen as the optimization goal. The Pythagorean line graph (PH) method is another impressive geometric method to solve the collaborative path planning problem [33,34,35]. Dynamic constraints are described using geometric differential properties such as curvature and torsion. Although these pure geometric methods have advantages for computational efficiency, they lack consideration of the dynamic characteristics of UAVs. In task allocation when designing numerous UAV trajectories, a practical task planning approach can effectively filter results by applying a pure geometric method. However, this approach may make some tracks difficult to fly. The trajectory generation method based on Bezier curves adopts a combination of geometric and dynamic methods, meaning it can be called a geometric–dynamic method [36]. Further, Petropoulos and Longuski introduced the shape function concept [37], and Jolly et al. used Bezier curves for robot path planning [38]. Bezier curves have been employed to link lines or circles to obtain better curvature smoothness and continuity for UAV trajectories [39,40]. Another similar approach is the B-spline method [41]. In terms of spacecraft trajectory planning, inspired by the shape function concept proposed by Petropoulos et al., Huo et al. applied Bezier curves to the initial three-dimensional trajectory planning of solar sails and proposed a method of quickly generating the minimum-time three-dimensional trajectories of spacecraft propelled by solar sails using feedback control. The simulation results show that, considering the actual characteristics of the thrust vector, this method can design the transfer trajectories of spacecraft propelled by solar sails equipped with the feedback control devices with high accuracy in approximately 1% of the time required by the traditional Gaussian pseudospectral method [42,43,44]. Xiaoliang studied methods of continuous-curvature bounded path planning of fixed-wing UAVs based on Bezier curves [45], and Yu et al. conducted research on the fast generation of the trajectories of multiple UAVs arriving simultaneously based on spatiotemporal Bezier curves [46]. In this study, we investigated the attitude–path coupling dynamics and optimal control of EAD-UAVs via the Bezier shaping approach (BSA).
The rest of this paper is organized as follows. Section 2 presents the EAD-UAV attitude–path coupling dynamics based on the one-dimensional dynamic model of an EAD thruster. Section 3 describes the implementation of the BSA and IBSA to realize rapid trajectory optimization considering the coupling dynamics constraints in order to realize the EAD-UAV flight control for specified targets. Section 4 presents the numerical results within an optimal framework using the proposed approach. Finally, Section 5 summarizes the conclusions.

2. EAD-UAV Attitude–Path Coupling Dynamics

2.1. Reference Frames

Describing different variables in different coordinate systems is necessary to study EAD-UAV flight dynamics more conveniently and to make the expressions of the dynamics equations clearer. Therefore, the coordinate systems used and the conversion relationships between these coordinate systems must be clarified and can be described as follows.
(1)
Ground coordinate system O G x G y G z G
The ground coordinate system, O G x G y G z G , is fixed to the ground, and the original point, O G , is selected as a point joined with the ground according to the application. The y G -axis points perpendicularly upward from the ground, the x G -axis points to the initial course of the EAD-UAV in the horizontal plane, and the z G -axis direction is determined by the right-hand rule. By moving the origin, O G , to coincide with the centroid O of the EAD-UAV, the ground coordinate system, O x g y g z g , which is convenient for coordinate system transformation, can be obtained.
(2)
Body coordinate system O x t y t z t
The body coordinate system, O x t y t z t , is fixed to the EAD-UAV where the origin O is located at its centroid, and the x t -axis points to the head of the EAD-UAV along the direction of the body axis. The y t -axis is located in the longitudinal symmetry plane of the EAD-UAV body, perpendicular to the x t -axis and pointing upwards. The z t -axis is perpendicular to the O x t y t plane, and the direction is determined by the right-hand rule.
(3)
Speed coordinate system O x V y V z V
The origin O is defined at the center of mass of the EAD-UAV. The x V -axis points in the direction of V U (the speed of the drone relative to the air). The y V -axis is located in the plane of longitudinal symmetry of the fuselage, perpendicular to the x V -axis and pointing upward. The z V -axis is perpendicular to the O x V y V plane, and the direction is determined by the right-hand rule.
(4)
Track coordinate system O x h y h z h
The origin O is defined at the center of mass of the EAD-UAV. The x h -axis points in the direction of V (the speed of the UAV relative to the ground). The y h -axis is perpendicular to the x h -axis, located in the vertical plane of V , pointing upward. The z h -axis is perpendicular to the O x h y h plane, and the direction is determined by the right-hand rule.
The transformation matrix, B g t , from O x g y g z g to O x t y t z t and the transformation matrix, B t g , from O x t y t z t to O x g y g z g are given by
B g t = B t g T = 1 0 0 0 cos γ sin γ 0 sin γ cos γ cos ϑ sin ϑ 0 sin ϑ cos ϑ 0 0 0 1 cos ψ 0 sin ψ 0 1 0 sin ψ 0 cos ψ ,
where the pitch angle, ϑ 90 , 90 , is the angle between the x t -axis and the O x g y g plane, the yaw angle, ψ 180 , 180 , is the angle between the projection of the x t -axis on the O x g y g plane and the x g -axis, and the roll angle, γ 180 , 180 , is the angle between the vertical plane where the x t axis is located and the O x t y t plane as shown in Figure 3.
The transformation matrix, B g h , from O x g y g z g to O x h y h z h and the transformation matrix, B h g , from O x h y h z h to O x g y g z g are given by
B g h = B h g T = cos θ sin θ 0 sin θ cos θ 0 0 0 1 cos ψ V 0 sin ψ V 0 1 0 sin ψ V 0 cos ψ V ,
where the track deflection angle, ψ V 180 , 180 , is the angle between the projection line of the x h -axis on the O x g y g plane and the x g -axis and the track inclination angle, θ 90 , 90 , is the angle between the x h -axis and the O x g y g plane as shown in Figure 4.
The transformation matrix, B V t , from O x V y V z V to O x t y t z t and the transformation matrix, B t V , from O x t y t z t to O x V y V z V are given by
B V t = B t V T = cos α sin α 0 sin α cos α 0 0 0 1 cos β 0 sin β 0 1 0 sin β 0 cos β ,
where the angle of attack, α 90 , 90 , is the angle between the projection of the x V -axis on the O x t y t plane and the x t -axis and the sideslip angle, β 180 , 180 , is the angle between the x V -axis and the O x t y t plane as shown in Figure 5.

2.2. Forces and Torques Acting on EAD-UAVs

The mechanical environments that EAD-UAVs face during flight are very complex. As a compromise between modeling accuracy and computational cost, only the following main forces and torques were considered in this study.
(1)
Aerodynamic forces
The directions of the aerodynamic forces acting on the EAD-UAV are described in O x V y V z V , where the aerodynamic forces can be decomposed into three components of drag (X), lift (Y), and lateral force (Z) along the three coordinate axes of O x V y V z V . The main factors influencing the aerodynamic forces are the dynamic pressure of the incoming flow, q V , and the characteristic area of the EAD-UAV, S . The relationships between the aerodynamic forces and influencing factors are given by
X = C X q V S Y = C Y q V S Z = C Z q V S q V = 1 2 ρ V U 2 ,
where C X , C Y , and C Z are the dimensionless drag coefficient, lift coefficient, and lateral force coefficient, respectively, and ρ is the air density.
(2)
Aerodynamic torques
The aerodynamic torques acting on an EAD-UAV affect its flight attitude and can be decomposed into three components: roll torque ( M x ), yaw torque ( M y ), and pitch torque ( M z ) relative to the three coordinate axes of O x t y t z t . In addition to the dynamic pressure and characteristic area, the characteristic lengths of the EAD-UAV are the main factors influencing the aerodynamic torques. The relationships between the aerodynamic torques and influencing factors are given by
M x = C m x q V S l c M y = C m y q V S l c M z = C m z q V S l z ,
where l c is the lateral characteristic length, l z is the longitudinal characteristic length, and C m x , C m y , and C m z are the dimensionless roll, yaw, and pitch moment coefficients, respectively.
(3)
EAD forces and torques
Unlike traditional fixed-wing UAVs, EAD-UAVs have no moving parts, which means they have lower mechanical fatigue, longer service lives, and lower flight noise. Therefore, the control torques cannot be obtained through the deflection of the rudder surface. In order to control the attitude of an EAD-UAV, it is necessary to control the thrust distribution by adjusting the voltage distribution in different areas of the six sets of thrusters installed orthogonally, so as to generate control torque and realize attitude control of the EAD-UAV. The x t -directional thrust, P x , and pitch torque, M z EAD , are generated by O z t x t plane-symmetric EAD thrusters 1 and 2. The z t -directional thrust, P z , and yaw torque, M y EAD , are generated by O y t z t plane-symmetric EAD thrusters 3 and 4. Finally, the y t -directional thrust, P y , and roll torque, M x EAD , are generated by O x t y t plane-symmetric EAD thrusters 5 and 6. The EAD thrusts and steering torques can be expressed in the form shown in Equation (6) using the thrusts of the six sets of EAD thrusters:
P x = F 1 + F 2 P y = F 5 + F 6 P z = F 3 + F 4 M x EAD = F 5 F 6 l 3 M y EAD = F 3 F 4 l 2 M z EAD = F 1 F 2 l 1 ,
where l 1 is the distance between the projection of the thrust center of thruster 1 (or thruster 2) on the O y t z t plane and the center of mass, l 2 is the distance between the projection of the thrust center of thruster 3 (or thruster 4) on the O x t y t plane and the centroid, and l 3 is the distance between the projection of the thrust center of thruster 5 (or thruster 6) on the O z t x t plane and the center of mass, as shown in Figure 2.

2.3. Dynamic Equations of EAD-UAVs

(1)
Dynamic equation of the motion of the center of mass of an EAD-UAV
In order to analyze the change law of an EAD-UAV track conveniently, the dynamic vector equation of the motion of the centroid can be projected into the track coordinate system as
m V ˙ V θ ˙ V ψ ˙ V cos θ = B g h B t g P x P y P z + B V t X Y Z + B g h 0 m g 0 ,
(2)
Dynamic equation of an EAD-UAV rotating around the centroid
The dynamic vector equation of an EAD-UAV rotating around the centroid can be projected into the body coordinate system as
J x ω ˙ x J x y ω ˙ y + J z J y ω z ω y + J x y ω x ω z J y ω ˙ y J x y ω ˙ x + J x J z ω x ω z J x y ω z ω y J z ω ˙ z + J y J x ω x ω y + J x y ω y 2 ω x 2 = M x + M x EAD M y + M y EAD M z + M z EAD ,
where ω x , ω y , and ω z are the angular velocities of the EAD-UAV rotating around each axis of O x t y t z t ; J x , J y , and J z are the moments of inertia of the EAD-UAV relative to each axis of the body coordinate system; and J x y is the inertial product of the EAD-UAV with respect to the x t - and y t -axes. Because the EAD-UAV is symmetrical on the longitudinal symmetry plane, O x t y t , and there is no symmetry deviation, J y z = J z x = 0 .
(3)
Kinematic equation of the motion of the center of mass of an EAD-UAV
According to the relationship between flight speed and position, the kinematic equation of the center of mass of an EAD-UAV in O G x G y G z G is given by
x ˙ y ˙ z ˙ = B h g V 0 0 = V cos θ cos ψ V V sin θ V cos θ sin ψ V .
(4)
Kinematic equation of an EAD-UAV rotating about its center of mass
In order to obtain the attitude information of an EAD-UAV in the air, it is necessary to establish a kinematic equation describing an EAD-UAV relative to O x g y g z g ; its form is
ϑ ˙ ψ ˙ γ ˙ = ω y sin γ + ω z cos γ 1 cos ϑ ω y cos γ ω z sin γ ω x tan ϑ ω y cos γ ω z sin γ .
(5)
Supplementary equations
When flying an EAD-UAV in non-calm atmosphere conditions, the effect of the wind speed on the flight must be considered. The relationship among the speed of the EAD-UAV relative to the ground, V , speed of the relative airflow, V U , and wind speed, V W , is
V = V U + V W ,
where V U can be expressed as shown in Equation (12); V U x t , V U y t , and V U z t are the V U -components of the three coordinate axes in O x t y t z t ; and V W x , V W y , and V Wz are the V W -components of the three coordinate axes in O G x G y G z G :
V U x t V U y t V U y t = B g t B h g V 0 0 V W x V W y V Wz .
The angle of attack α and sideslip angle β of an EAD-UAV can be expressed as
tan α = V U y t / V U x t sin β = V U z t / V U V U = V U x t 2 + V U y t 2 + V U z t 2 ,
(6)
Relationship between the thrust and voltage of an EAD thruster
In the one-dimensional model of an EAD thruster cited from [15], its thrust is given by
F = I d μ = C U U U 0 d μ ,
where I is the current of the EAD thruster; C is a constant related to the EAD thruster structure and ion mobility, μ ; U is the applied voltage; U 0 is the initial voltage of the EAD thruster when corona discharge occurs; and d is the distance between the emitter electrode and collector electrode of the EAD thruster electrode pair. C can be expressed as Equation (15), cited from [14],
C = C 0 l μ ε 0 d 2 ,
where l is the electrode length of the EAD thruster and ε 0 is the vacuum permittivity. The dimensionless constant, C 0 , is estimated based on the relevant data in [9].
U 0 can be determined applying Equation (15), cited from [47], by approximating the situation in this paper to the form of parallel wires,
U 0 = E p r c ln d r c .
Here, E p is the intensity of the corona inception electric field at which the corona discharge phenomenon occurs, which can be calculated by utilizing Equation (17), cited from [48], and r c is the radius of the emitter electrode in centimeters,
E p = E 0 δ ε 1 + 0.308 δ r c .
Here, E 0 is the electric field strength corresponding to air breakdown (for air, E 0 = 31 kV/cm) δ is the relative atmospheric density, and ε is the smoothness of the electrode surface; in this paper, ε = 1 .
(7)
Coupling dynamic equation of an EAD-UAV
The dynamic equation of an EAD-UAV can be expressed in the form of x ˙ = f ( x , u ) , where x = x , y , z , ϑ , ψ , γ , ω x , ω y , ω z , V , θ , ψ V T is the state variable vector and u = U 1 , U 2 , U 3 , U 4 , U 5 , U 6 T is the control variable vector.
According to the above equations, the position–attitude coupling dynamic equation of the EAD-UAV can be obtained, as shown in Equation (18). The coupled control of the position and attitude of the EAD-UAV can be achieved by adjusting the voltages of the six EAD thrusters.
x ˙ y ˙ z ˙ ϑ ˙ ψ ˙ γ ˙ ω ˙ x ω ˙ y ω ˙ z V ˙ θ ˙ ψ ˙ V = V cos θ cos ψ V V sin θ V cos θ sin ψ V ω y sin γ + ω z cos γ 1 cos ϑ ω y cos γ ω z sin γ ω x tan ϑ ω y cos γ ω z sin γ J x y 2 ω y ω z + J x + J y J z ω x ω z F 3 u F 4 u l 2 M y J x y J y J z ω y ω z + F 5 u F 6 u l 3 + M x J y J x J y + J x y 2 J x y 2 ω x ω z + J x J y + J z ω y ω z F 5 u F 6 u l 3 M x J x y J x + J z ω x ω z + F 3 u F 4 u l 2 + M y J x J x J y + J x y 2 J x y ω x 2 + J x J y ω x ω y J x y ω y 2 + F 1 u F 2 u l 1 + M z J z A 4 A 5 A 6 cos ψ A 7 cos θ A 3 sin ϑ A 1 cos ϑ cos γ + A 2 cos ϑ sin γ + G sin θ m A 4 + A 5 + A 6 cos ψ V + A 7 sin θ A 3 sin ϑ A 1 cos ϑ cos γ + A 2 cos ϑ sin γ + G cos θ m V A 2 cos γ A 1 sin γ cos ψ + A 8 sin ψ cos ψ V A 8 cos ψ + A 2 cos γ + sin γ A 1 sin ψ sin ψ V m V cos θ ,
where
A 1 = Z sin β + X cos β sin α + Y cos α + F 5 u + F 6 u A 2 = Z cos β X sin β + F 3 u + F 4 u A 3 = Z sin β cos α + X cos α cos β Y sin α F 1 u F 2 u A 4 = A 1 cos γ + A 2 sin γ cos ψ cos ψ V + sin ψ sin ψ V sin ϑ A 5 = A 2 cos ψ sin ψ V sin ψ cos ψ V cos γ A 6 = A 3 cos ϑ cos ψ A 1 sin γ sin ψ A 7 = A 1 sin γ cos ψ + A 3 sin ψ cos ϑ sin ψ V A 8 = A 1 cos γ + A 2 sin γ sin ϑ A 3 cos ϑ F i ( u ) = C U i U i U 0 d μ , i = 1 , 2 , 3 , 4 , 5 , 6

3. Trajectory Optimization Using the Integrative Bezier Shaping Approach (IBSA)

3.1. Optimization Problem

This paper describes the trajectory planning problem of multi-point exploration using an EAD-UAV. The EAD-UAV leaves the starting point of the flight trajectory at moment T 0 and reaches the first target point at moment T 1 . Then, the EAD-UAV arrives at the k-th objective at moment T k , where the superscript k = [1, N] is the parameter and corresponds to the k-th segment of the trajectory (that is, the trajectory between the (k–1)-th objective and the k-th objective). The optimal EAD-UAV control law can be obtained by taking the minimum total flight time of the flight trajectory of an EAD-UAV flying over the N targets as the optimization objective, that is, the voltage distribution of the six EAD thrusters at each moment, u = U 1 U 2 U 3 U 4 U 5 U 6 T . Therefore, the objective function of the optimization problem can be expressed as
J T = k = 1 N Δ T k ,
where Δ T k = T k T k 1 is the flight time from the (k–1)-th objective to the k-th objective. The flight times Δ T k , k = [ 1 , N ] are chosen as the optimization variables in the process of solving.
EAD-UAV has high voltage requirements, which places high demands on batteries, so minimum energy consumption is also a reasonable optimization objective. The objective function for minimum energy consumption can be expressed as:
J E n e r g y = k = 1 N j = 1 m 1 i = 1 6 U i j I i j Δ T k j ,
where Δ T k j is the time interval between the j-th point and the (j + 1)-th point after discretizing the time between T k 1 and T k into m points, while P p o w e r = U i j I i j is the power of the i-th thruster at the j-th point. P a v e r a g e = J E n e r g y Δ T k is the average power of the k-th trajectory. The average power can be used to evaluate the battery performance requirements of EAD-UAV.
For the multi-point continuous optimization problem of EAD-UAV trajectories, besides the boundary constraints at moments T 0 and T N , the boundary constraints of the k-th trajectory at moment T k ,   k = [ 1 , N 1 ] should be considered. The position, velocity, track angles, and attitude angles of the EAD-UAV in the ground coordinate system at T 0 and T N and the triaxial angular velocity in the body coordinate system at T 0 and T N are known quantities; the boundary constraints at moments T 0 and T N are given by
x T 0 = x T 0 y T 0 z T 0 V T 0 θ T 0 ψ V T 0    ϑ T 0 ψ T 0 γ T 0 ω x T 0 ω y T 0 ω z T 0 T x T N = x T N y T N z T N V T N θ T N ψ V T N    ϑ T N ψ T N γ T N ω x T N ω y T N ω z T N T ,
where x is the state vector of the EAD-UAV.
At moment T k ,   k = [ 1 , N 1 ] , the positions are known quantities as boundary constraints. Therefore, the boundary constraints are given by
x T k = x k T k y T k = y k T k z T k = z k T k ,
where x k T k , y k T k , and z k T k are the three-axis positions of the k-th objective at time T k , k = [ 1 , N 1 ] .
Furthermore, although there are no boundary constraints for the velocity, acceleration, attitude angle, or angular velocity at the intermediate points, these quantities also need to be continuous for the EAD-UAV trajectory to be continuous.

3.2. Bezier State Approximation

Among the state variables in the EAD-UAV dynamics equation, ω x , ω y , ω z , V , θ , ψ V can be expressed as follows by using x , y , z , ϑ , ψ , γ and their first derivatives:
ω x ω y ω z V θ ψ V = γ ˙ + ψ ˙ cos ϑ tan ϑ ϑ ˙ sin γ + ψ ˙ cos γ cos ϑ ϑ ˙ cos γ ψ ˙ sin γ cos ϑ x ˙ 2 + y ˙ 2 + z ˙ 2 arcsin y ˙ V arctan z ˙ x ˙ .
Therefore, all state variables and their first derivatives can be represented by x , y , z , ϑ , ψ , γ and their first and second derivatives.
According to the EAD-UAV dynamic equation, the thrusts F 1 F 6 of the six EAD thrusters can be determined and expressed as follows by using x , y , z , ϑ , ψ , γ and their first and second derivatives:
F 1 = m l 1 B 1 B 2 cos ϑ + B 3 + B 4 + B 5 2 l 1 F 2 = m l 1 B 1 B 2 cos ϑ + B 3 + B 4 B 5 2 l 1 F 3 = m l 2 Γ 1 Γ 2 cos θ Γ 3 + Γ 4 + Γ 5 2 l 2 F 4 = m l 2 Γ 1 Γ 2 cos θ Γ 3 + Γ 4 Γ 5 2 l 2 F 5 = l 3 Δ 1 Δ 2 cos γ Δ 3 + Δ 4 + Δ 5 2 l 3 F 6 = l 3 Δ 1 Δ 2 cos γ Δ 3 + Δ 4 Δ 5 2 l 3 ,
where
B 1 = ψ V sin ψ V + V ˙ cos ψ cos ψ V + sin ψ V ψ ˙ V cos ψ V + V ˙ sin ψ cos θ B 2 = θ ˙ sin θ V cos ψ cos ψ V + sin ψ sin ψ V B 3 = l 1 V ϑ m cos θ sin ϑ + l 1 V ˙ sin θ m + G sin ϑ B 4 = l 1 Z sin β + X cos β cos α Y l 1 sin α B 5 = J x y ω x 2 ω y J x J y ω x + J x y ω y 2 + J z ω ˙ z M z Γ 1 = ψ ˙ V sin ψ V + V ˙ cos ψ cos ψ V + sin ψ V ψ ˙ V cos ψ V + V ˙ sin ψ sin ϑ ϑ ˙ cos ϑ V sin γ Γ 2 = ψ ˙ V cos ψ V V ˙ sin ψ cos ψ V + sin ψ V ψ ˙ V sin ψ V + V ˙ cos ψ cos γ Γ 3 = θ ˙ sin θ V m cos ψ cos ψ V + sin ψ sin ψ V sin ϑ + cos ϑ V ˙ sin θ m + G l 2 sin γ Γ 4 = θ ˙ sin θ V m l 2 cos ψ sin ψ V cos ψ V sin ψ cos γ Z l 2 cos β + X l 2 sin β Γ 5 = ω y J x y + ω x J x J z ω z J x y ω ˙ x + J y ω ˙ y M y Δ 1 = m ψ ˙ V sin ψ V + V ˙ cos ψ cos ψ V + sin ψ V ψ ˙ V cos ψ V + V ˙ sin ψ sin ϑ θ ˙ cos ϑ V cos θ Δ 2 = θ ˙ sin θ V m cos ψ cos ψ V + sin ψ sin ψ V sin ϑ cos ϑ V ˙ sin θ m + G Δ 3 = m l 3 sin γ ψ ˙ V cos ψ V V ˙ sin ψ cos ψ V + sin ψ V ψ ˙ V sin ψ V + V ˙ cos ψ cos θ Δ 4 = l 3 θ ˙ sin θ V m cos ψ sin ψ V cos ψ V sin ψ sin γ l 3 Z sin β + X cos β sin α Y l 3 cos α Δ 5 = ω x J x y ω y J y J z ω z J x y ω ˙ y + J x ω ˙ x M x
It can be seen from Equation (14) that after corona discharge occurs (that is, when the applied voltage is greater than the corona initiation voltage), the thrust of the EAD thruster increases monotonically with increasing voltage, so the applied voltage, U , can be given by Equation (24) in the form of a function of the EAD thrust, F . Therefore, the control variables of the dynamic equation of an EAD-UAV can be expressed by using x , y , z , ϑ , ψ , γ and their first and second derivatives.
U = C U 0 + C C U 0 2 + 4 F 2 C ,
Here, C is given by
C = C d μ .
Therefore, the whole EAD-UAV dynamic equation can be expressed by using x , y , z , ϑ , ψ , γ and their first and second derivatives.
Concerning the k-th trajectory, the approximations of the six-degree-of-freedom position can be expanded in the time domain by applying a Bezier curve function:
x ( k ) ( τ ) = j = 0 n x ( k ) B x , j ( k ) ( τ ) P x , j ( k ) y ( k ) ( τ ) = j = 0 n y ( k ) B y , j ( k ) ( τ ) P y , j ( k ) z ( k ) ( τ ) = j = 0 n z ( k ) B z , j ( k ) ( τ ) P z , j ( k ) ϑ ( k ) ( τ ) = j = 0 n ϑ ( k ) B ϑ , j ( k ) ( τ ) P ϑ , j ( k ) ψ ( k ) ( τ ) = j = 0 n ψ ( k ) B ψ , j ( k ) ( τ ) P ψ , j ( k ) γ ( k ) ( τ ) = j = 0 n γ ( k ) B γ , j ( k ) ( τ ) P γ , j ( k ) ,
where 0 τ = t / Δ T ( k ) 1 is the self-defined scaled time parameter; t and Δ T ( k ) are the actual time and total flight time duration in the k-th segment of trajectory, respectively; n x ( k ) , n y ( k ) , n z ( k ) , n ϑ ( k ) , n ψ ( k ) , and n γ ( k ) are the orders of the Bezier curve function for each of the six degrees; and P x , j ( k ) ( j [ 0 , n x ( k ) ] ) , P y , j ( k ) ( j [ 0 , n y ( k ) ] ) , P z , j ( k ) ( j [ 0 , n z ( k ) ] ) , P ϑ , j ( k ) ( j [ 0 , n ϑ ( k ) ] ) , P ψ , j ( k ) ( j [ 0 , n ψ ( k ) ] ) , and P γ , j ( k ) ( j [ 0 , n ψ ( k ) ] ) are the Bezier coefficients. B x , j ( k ) ( τ ) , B y , j ( k ) ( τ ) , B z , j ( k ) ( τ ) , B ϑ , j ( k ) ( τ ) , B ψ , j ( k ) ( τ ) , and B γ , j ( k ) ( τ ) are the Bezier basis functions for the coordinate approximations and are given by
B x , j ( k ) ( τ ) = n x ( k ) ! j ! n x ( k ) j ! τ j ( 1 τ ) n x ( k ) j j [ 0 , n x ( k ) ] B y , j ( k ) ( τ ) = n y ( k ) ! j ! n y ( k ) j ! τ j ( 1 τ ) n y ( k ) j j [ 0 , n y ( k ) ] B z , j ( k ) ( τ ) = n z ( k ) ! j ! n z ( k ) j ! τ j ( 1 τ ) n z ( k ) j j [ 0 , n z ( k ) ] B ϑ , j ( k ) ( τ ) = n ϑ ( k ) ! j ! n ϑ ( k ) j ! τ j ( 1 τ ) n ϑ ( k ) j j [ 0 , n ϑ ( k ) ] B ψ , j ( k ) ( τ ) = n ψ ( k ) ! j ! n ψ ( k ) j ! τ j ( 1 τ ) n ψ ( k ) j j [ 0 , n ψ ( k ) ] B γ , j ( k ) ( τ ) = n γ ( k ) ! j ! n γ ( k ) j ! τ j ( 1 τ ) n γ ( k ) j j [ 0 , n γ ( k ) ] ,
In order to avoid repetition, individually, the procedure of an approximation for coordinate x is presented in this paper, and the other five degrees can be handled in an analogical manner. According to Equation (27), the first and second τ-derivatives of the degree x can be written as
x ( k ) ( τ ) = j = 0 n x ( k ) B x , j ( k ) ( τ ) P x , j ( k ) x ( k ) ( τ ) = j = 0 n x ( k ) B x , j ( k ) ( τ ) P x , j ( k ) ,
where the superscript ′ represents the derivative with respect to the parameter τ and B x , j ( k ) ( τ ) and B x , j ( k ) ( τ ) are the first and second τ-derivatives of the Bezier basis function for the degree x approximation in the k-th interplanetary trajectory, respectively. According to Equation (28), one can obtain
B x , j ( k ) ( τ ) = n x ( k ) ( 1 τ ) n x ( k ) 1 j = 0 n x ( k ) ! ( j 1 ) ! n x ( k ) j ! τ j 1 ( 1 τ ) n x ( k ) j n x ( k ) ! j ! n x ( k ) j 1 ! τ j ( 1 τ ) n x ( k ) j 1 j [ 1 , n x ( k ) 1 ] n x ( k ) τ n x ( k ) 1 j = n x ( k ) B x , j ( k ) ( τ ) = n x ( k ) n x ( k ) 1 ( 1 τ ) n x ( k ) 2 j = 0 n x ( k ) n x ( k ) 1 n x ( k ) 2 τ ( 1 τ ) n x ( k ) 3 2 n x ( k ) n x ( k ) 1 ( 1 τ ) n x ( k ) 2 j = 1 n x ( k ) ! ( j 2 ) ! n x ( k ) j ! τ j 2 ( 1 τ ) n x ( k ) j 2 n x ( k ) ! ( j 1 ) ! n x ( k ) j 1 ! τ j 1 ( 1 τ ) n x ( k ) j 1 + n x ( k ) ! j ! n x ( k ) j 2 ! τ j ( 1 τ ) n x ( k ) j 2 j [ 2 , n x ( k ) 2 ] n x ( k ) n x ( k ) 1 n x ( k ) 2 τ n x ( k ) 3 ( 1 τ ) 2 n x ( k ) n x ( k ) 1 τ n x ( k ) 2 j = n x ( k ) 1 n x ( k ) n x ( k ) 1 τ n x ( k ) 2 j = n x ( k ) ,
By substituting τ = 0 and τ = 1 into Equations (29) and (30), interesting characteristics of the Bezier basis function at the boundary can be obtained, as follows:
B x , j ( k ) ( τ = 0 ) = 1 j = 0 0 j [ 1 , n x ( k ) ] B x , j ( k ) ( τ = 1 ) = 0 j [ 0 , n x ( k ) 1 ] 1 j = n x ( k ) B x , j ( k ) ( τ = 0 ) = n x ( k ) j = 0 n x ( k ) j = 1 0 j [ 2 , n x ( k ) ] B x , j ( k ) ( τ = 1 ) = 0 j [ 0 , n x ( k ) 2 ] n x ( k ) j = n x ( k ) 1 n x ( k ) j = n x ( k ) B x , j ( k ) ( τ = 0 ) = n x ( k ) n x ( k ) 1 j = 0 2 n x ( k ) n x ( k ) 1 j = 1 n x ( k ) n x ( k ) 1 j = 2 0 j [ 3 , n x ( k ) ] B x , j ( k ) ( τ = 1 ) = 0 j [ 0 , n x ( k ) 3 ] n x ( k ) n x ( k ) 1 j = n x ( k ) 2 2 n x ( k ) n x ( k ) 1 j = n x ( k ) 1 n x ( k ) n x ( k ) 1 j = n x ( k ) ,
According to Equations (27), (29) and (31), the relationships between the boundary constraints of the k-th segment of the trajectory and the Bezier coefficient can be found, as shown in Equation (32):
x T ( k 1 ) = x ( k ) ( τ = 0 ) = P x , 0 ( k ) x T ( k ) = x ( k ) ( τ = 1 ) = P x , n x ( k ) ( k ) Δ T ( k ) x ˙ T ( k 1 ) = x ( k ) ( τ = 0 ) = n x ( k ) P x , 1 ( k ) P x , 0 ( k ) Δ T ( k ) x ˙ T ( k ) = x ( k ) ( τ = 1 ) = n x ( k ) P x , n x ( k ) ( k ) P x , n x ( k ) 1 ( k ) Δ T ( k ) 2 x ¨ T ( k 1 ) = x ( k ) ( τ = 0 ) = n x ( k ) n x ( k ) 1 P x , 0 ( k ) + P x , 2 ( k ) 2 P x , 1 ( k ) Δ T ( k ) 2 x ¨ T ( k ) = x ( k ) ( τ = 1 ) = n x ( k ) n x ( k ) 1 P x , n x ( k ) ( k ) + P x , n x ( k ) 2 ( k ) 2 P x , n x ( k ) 1 ( k ) ,
where the superscript · represents the derivative of the actual time, t, from which one can obtain the Bezier representation of x and its first and second derivatives at the start time T ( k 1 ) and end time T ( k ) of the k-th segment of the trajectory.
Hence, it was ensured that the approximate EAD-UAV flight trajectory naturally met the boundary constraints in aspects of position, velocity, and acceleration at the boundary time ( T ( k 1 ) and T ( k ) ) of the k-th segment by determining Bezier coefficients P x , 0 ( k ) , P x , 1 ( k ) , P x , 2 ( k ) , P x , n x ( k ) 2 ( k ) , P x , n x ( k ) 1 ( k ) , and P x , n x ( k ) ( k ) as
P x , 0 ( k ) = x T ( k 1 ) P x , 1 ( k ) = x T ( k 1 ) + Δ T ( k ) x ˙ T ( k 1 ) / n x ( k ) P x , 2 ( k ) = x T ( k 1 ) + 2 Δ T ( k ) x ˙ T ( k 1 ) / n x ( k ) + Δ T ( k ) 2 x ¨ T ( k 1 ) / n x ( k ) n x ( k ) 1 P x , n x ( k ) 2 ( k ) = x T ( k ) 2 Δ T ( k ) x ˙ T ( k ) / n x ( k ) + Δ T ( k ) 2 x ¨ T ( k ) / n x ( k ) n x ( k ) 1 P x , n x ( k ) 1 ( k ) = x T ( k ) Δ T ( k ) x ˙ T ( k ) / n x ( k ) P x , n x ( k ) ( k ) = x T ( k ) ,
This characteristic is fairly beneficial for dealing with fast trajectory optimization; instead of concerning boundary constraints, only the dynamic constraints and optimization performance index must be considered during the optimization procedure. In order to deal with the acceleration continuity of an EAD-UAV in multi-target trajectory optimization, the relationship between the second derivative of the coordinate approximations ( x ¨ T ( k 1 ) , x ¨ T ( k ) ) and the coefficients ( P x , 2 ( k ) , P x , n x ( k ) 2 ( k ) ) is determined, which is different from the single-target trajectory optimization problem.

3.3. Nonlinear Programming Problem (NLP)

The essential problem is to obtain a practical trajectory within the capabilities of the EAD thrusters in all segments of EAD-UAV flight. Thus, it is essential to assess the motion of each trajectory at selected discretized points. Legendre–Gauss (LG) distribution, defined as the roots of the m ( k ) th-degree Legendre polynomial, was selected as the discrete mode of discretized points in this paper and is given by
τ 1 = 0 < τ 2 < < τ m ( k ) 1 < τ m ( k ) = 1 ,
The discrete six-degree-of-freedom position and attitude coordinates can be expressed in the form of matrix products. For example, x can be expressed as
[ x ( k ) ] m ( k ) × 1 = [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) [ P x ( k ) ] ( n x ( k ) + 1 ) × 1 [ x ( k ) ] m ( k ) × 1 = [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) [ P x ( k ) ] ( n x ( k ) + 1 ) × 1 [ x ( k ) ] m ( k ) × 1 = [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) [ P x ( k ) ] ( n x ( k ) + 1 ) × 1 ,
where [ P x ( k ) ] ( n x ( k ) + 1 ) × 1 = [ P x , 0 ( k ) P x , 1 ( k ) P x , n x ( k ) 1 ( k ) P ρ , n x ( k ) ( k ) ] T is a column vector composed of known and unknown Bezier coefficients and P x , 0 ( k ) , P x , 1 ( k ) , P x , 2 ( k ) , P x , n x ( k ) 2 ( k ) , P x , n x ( k ) 1 ( k ) , P x , n x ( k ) ( k ) are determined by the boundary constraints. At the middle target point of the trajectory, the boundary condition is only the position of the target point. Considering the continuity of velocity and acceleration, P x , 0 ( k ) , P x , 1 ( k ) , P x , 2 ( k ) can be determined according to
P x , 0 ( k ) = P x , n x ( k 1 ) ( k 1 ) P x , 1 ( k ) = 1 + Δ T ( k ) Δ T ( k 1 ) P x , n x ( k 1 ) ( k 1 ) Δ T ( k ) Δ T ( k 1 ) P x , n x ( k 1 ) 1 ( k 1 ) P x , 2 ( k ) = 1 + Δ T ( k ) Δ T ( k 1 ) 2 P x , n x ( k 1 ) ( k 1 ) 2 Δ T ( k ) Δ T ( k 1 ) 1 + Δ T ( k ) Δ T ( k 1 ) P x , n x ( k 1 ) 1 ( k 1 ) + Δ T ( k ) Δ T ( k 1 ) 2 P x , n x ( k 1 ) 2 ( k 1 ) ,
By substituting discretized points, [ τ ] m ( k ) × 1 , into the coordinate approximations and their τ-derivatives, matrices [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) , [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) , and [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) can be expressed as follows:
[ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) = B x , 0 ( k ) ( τ 1 ) B x , n x ( k ) ( k ) ( τ 1 ) B x , 0 ( k ) ( τ m ( k ) ) B x , n x ( k ) ( k ) ( τ m ( k ) ) [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) = B x , 0 ( k ) ( τ 1 ) B x , n x ( k ) ( k ) ( τ 1 ) B x , 0 ( k ) ( τ m ( k ) ) B x , n x ( k ) ( k ) ( τ m ( k ) ) [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) = B x , 0 ( k ) ( τ 1 ) B x , n x ( k ) ( k ) ( τ 1 ) B x , 0 ( k ) ( τ m ( k ) ) B x , n x ( k ) ( k ) ( τ m ( k ) ) ,
It is worth noting that when the number of discretized points, m ( k ) , and the order of the Bezier curve, n x ( k ) , are determined, [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) , [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) , and [ B x ( k ) ] m ( k ) × ( n x ( k ) + 1 ) will be determined and become constant matrices. Thus, without the need for repeating calculation of the basis function matrices, the computational efficiency of the trajectory optimization greatly improves.
By using the six-degree-of-freedom coordinates in matrix form and their first and second τ -derivatives, the voltages of the six groups of plasma thrusters can be expressed in the form of the following discretization matrix:
[ U i ( k ) ] m ( k ) × 1 = U i ( k ) [ x ( k ) ] m ( k ) × 1 , [ y ( k ) ] m ( k ) × 1 , [ z ( k ) ] m ( k ) × 1 , [ ϑ ( k ) ] m ( k ) × 1 , [ ψ ( k ) ] m ( k ) × 1 , [ γ ( k ) ] m ( k ) × 1 , [ x ( k ) ] m ( k ) × 1 , [ y ( k ) ] m ( k ) × 1 , [ z ( k ) ] m ( k ) × 1 , [ ϑ ( k ) ] m ( k ) × 1 [ ψ ( k ) ] m ( k ) × 1 , [ γ ( k ) ] m ( k ) × 1 , [ x ( k ) ] m ( k ) × 1 , [ y ( k ) ] m ( k ) × 1 , [ z ( k ) ] m ( k ) × 1 , [ ϑ ( k ) ] m ( k ) × 1 , [ ψ ( k ) ] m ( k ) × 1 , [ γ ( k ) ] m ( k ) × 1 , ( i = 1 , 2 , 3 , 4 , 5 , 6 ) ,
When the EAD thruster is working, the voltage, U , must be greater than the corona voltage, U 0 . However, if the voltage is higher than the breakdown voltage, spark discharge will occur and the thruster will fail. Thus, the voltage of the thruster needs to satisfy the constraint U 0 U i U max , ( i = 1 , 2 , 3 , 4 , 5 , 6 ) . It is assumed that the drag coefficient, C X , and lift coefficient, C Y , remain unchanged during the flight of the EAD-UAV, so the angle of attack and sideslip angle are limited to a relatively small range, between α [ 1 , 1 ] and β [ 1 , 1 ] . Combined with the objective function shown in Equation (22), the N-segment trajectory optimization problem can be converted into a nonlinear programming problem (NLP) expressed as:
min [ X x ( k ) ] , [ X y ( k ) ] , [ X z ( k ) ] , [ X ϑ ( k ) ] , [ X ψ ( k ) ] , [ X γ ( k ) ] , Δ T ( k ) k = 1 N Δ T ( k ) s . t .    U 0 U i U max , ( i = 1 , 2 , 3 , 4 , 5 , 6 ) α min α α max β min β β max ,
where [ X x ( k ) ] , [ X y ( k ) ] , [ X z ( k ) ] , [ X ϑ ( k ) ] , [ X ψ ( k ) ] , [ X γ ( k ) ] are the unknown parts of the Bezier coefficients P x ( k ) , P y ( k ) , P z ( k ) , P ϑ ( k ) , P ψ ( k ) , P γ ( k ) , which need to be optimized to satisfy the dynamic constraints and obtain the optimal parameters.
The NLP whose optimization objective is minimum energy consumption can be expressed as:
min [ X x ( k ) ] , [ X y ( k ) ] , [ X z ( k ) ] , [ X ϑ ( k ) ] , [ X ψ ( k ) ] , [ X γ ( k ) ] , Δ T ( k ) k = 1 N j = 1 m 1 i = 1 6 U i ( j ) I i ( j ) Δ T ( k ) ( j ) s . t . U 0 U i U max , ( i = 1 , 2 , 3 , 4 , 5 , 6 ) α min α α max β min β β max ,

4. Numerical Results

An EAD-UAV was applied in two flight scenarios in this study: single- and multi-target continuous optimal flight control. The position and attitude parameters of the EAD-UAV at the starting and target points are given in Table A2 and Table A3. In the single-target optimal flight scenario, two numerical simulations for time–optimal trajectory optimization were conducted and the optimized results achieved by the BSA were used as feasible initial value estimations of the Gaussian pseudospectral method (GPM). The BSA is cited from [49] and the GPM is cited from [50]. Optimization for optimal trajectory of energy consumption was conducted using BSA, the results were compared with those obtained in time–optimal trajectory optimization. In the multi-target optimal flight control scenario, to demonstrate the practicability of the IBSA, a series of numerical simulations were conducted, and the results were compared with those obtained using the traditional BSA, which converts the optimization problem of multi-target EAD-UAV flight trajectory into multiple independently solved NLP problems. In the trajectory optimization problems using the IBSA and BSA, the Bezier orders and the number of LG points were chosen as n x ( k ) = n y ( k ) = n z ( k ) = n ϑ ( k ) = n ψ ( k ) = n γ ( k ) = 6 and m ( k ) = 50 . The interior-point method was employed to solve the converted NLP problem because of the absence of equality constraints. In the trajectory optimization problems using the GPM, the SQP algorithm was applied to solve the NLP problem, and the numbers of LG points were chosen to be 70. The interior-point method and the SQP algorithm were implemented by using the fmincon MATLAB function. All numerical simulations were performed on a Ryzen R7-5800H CPU 3.5 GHz with Windows 11 and run on MATLAB R2020a.

4.1. Single-Target Optimal Flight Control

To verify the practicability of the EAD-UAV, we conducted two numerical simulations for the time optimal single-target flight control of the BSA and GPM under different maximum voltages, Umax. Optimal single-target trajectories of the EAD-UAV with Umax = 80 kV were designed by the BSA and GPM and are shown in Figure 6. The mass of the EAD-UAV was 2.6 kg, and when U0 = 7.7 kV and Umax = 80 kV, the maximum thrust of each EAD thruster, Fmax, was 14.4184 N. Except for the boundary constraints between the start and target, there were no constraints on the speed and attitude of the EAD-UAV, only on the Umax of the EAD thruster. It can be seen from Figure 7 that the attitude angles of the EAD-UAV are within reasonable ranges during flight. As can be seen from Figure 8, the angle of attack and sideslip angle of the EAD-UAV are within the limited range, which satisfies the assumption that the lift coefficient and drag coefficient are constant. Figure 9 and Figure 10 reveal that the thrust of the EAD thruster always satisfies the maximum voltage constraint during flight, indicating that the dynamic model of the EAD-UAV is reasonable and that the BSA can fully satisfy the dynamic constraints. It can be seen from Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 that the BSA and GPM optimization results are very similar. The flight time of the optimal trajectory based on the BSA is 192.5029 s, and the calculation time is 0.7824 s. Meanwhile, the flight time of the optimal trajectory obtained through GPM is 189.8874 s, and the calculation time is 49.3636 s. The difference between the BSA and GPM results is 1.38%, and the calculation time of the BSA is only 1.58% of that of GPM. Thus, the BSA has considerable advantages in trajectory planning. In the figures, the +/− of U and F represent the thrust direction, + indicates that the EAD thruster generates the right or upward thrust, and − indicates the opposite meaning.
Table A4 presents the flight times obtained from the numerical simulations of the trajectory optimization using the BSA and GPM when Umax was in the range of [50:2:80] kV. With increasing Umax, the flight time decreases. A larger Umax makes the EAD-UAV have a higher flight speed and greater maneuverability. This result is in agreement with the expectations. When the voltage is higher than 50 kV and the thrust is higher than 5.2735 N, the simulation can converge to get the optimal trajectory. This level of thrust is similar to that in [9], and we can optimize the EAD thruster to reach this level, which shows that it is feasible to use BSA to optimize the trajectory of an EAD-UAV. The average deviation between the flight times of the optimal trajectories obtained by the BSA and GPM is 1.14%, and the average calculation time of the BSA is only 1.95% of that of GPM, which shows that the BSA has obvious advantages in terms of calculation efficiency and that the calculation accuracy is not much different from that of GPM.
To verify the battery performance requirements of EAD-UAV, we conducted a numerical simulation to optimize the trajectory under different voltages for the optimal energy consumption using BSA. Optimal single-target trajectories of the EAD-UAV with Umax = 60 kV and Umax = 80 kV designed by the BSA are shown in Figure 11. It can be seen from Figure 12 that the attitude angles of the EAD-UAV are within reasonable ranges during flight. As can be seen from Figure 13, the angle of attack and sideslip angle of the EAD-UAV are within the limited range. Figure 14 and Figure 15 reveal that the thrust of the EAD thruster always satisfies the maximum voltage constraint during flight. The maximum thrust during flight is only 4.2 N, which is similar to that in [9]. This level of thrust is expected to be achieved by improving the EAD thruster. It can be seen from Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 that the trajectory and other indexes under different voltages are very close, indicating that the voltage required for flight under the optimal energy consumption is relatively low. As shown in Table A5, the flight time, energy consumption, and average power of the flight trajectory with the optimal energy consumption were basically the same when Umax was in the range of [60:2:80], indicating that applied voltage was much lower than Umax to reduce energy consumption, leading to relatively lower battery performance requirements. The UAV power in [9] is 600 W, and the average power of the EAD-UAV under the optimal energy consumption is 1.4 kW. Considering that the EAD-UAV has more thrusters, the higher power is reasonable. Battery performance can be optimized based on the UAV in [9] to achieve this level of power. In the case of optimal time, the thrusters tend to increase the power as much as possible to reduce the flight time. The average power was 4~8 times that in the case of the optimal energy consumption, reducing the flight time by 18–28%. When Umax was in the range of [60:2:80], the average power reached 5~11 kW, which is very demanding for the battery, meaning it may be difficult to find a qualified power supply. The energy consumption taking minimum energy consumption as the optimization objective is 106.6640 W·h, much less than that in the case of optimal time, which can be used as a reference for the selection of battery capacities of EAD-UAVs in the design of power supply systems. The energy consumption taking time as the optimization objective is 323.9982~589.3262 W·h, which is 3~5.5 times that in the case of the optimal energy consumption. Although higher energy consumption can achieve rapid transfer between targets, larger energy consumption requires greater battery capacity. This will lead to greater battery weight, which will pose a greater challenge to the design of the EAD-UAV, considering its low thrust. Therefore, balancing the energy consumption and flight time is very important when optimizing the flight trajectory of EAD-UAVs.

4.2. Multi-Target Continuous Optimal Flight Control

To demonstrate the feasibility of the EAD-UAV and computational efficiency of the IBSA in multi-target continuous flight trajectory optimization, we analyzed the continuous trajectory optimization problem involving an EAD-UAV passing through three targets from the start. The IBSA and BSA were used for trajectory optimization. The acceleration continuity problem at the target point of the IBSA was taken as the constraint in the overall trajectory planning. The BSA performed subsection optimization of the three-stage trajectory and added the acceleration continuity constraint at the trajectory connection. The optimal flight path obtained using the IBSA when the maximum voltage of the EAD thruster was 80 kV is shown in Figure 16. There is a smooth transition at the connection of the three-segment trajectory, and the flight trajectory is a smooth curve. It can be seen from Figure 17 and Figure 18 that the velocity and acceleration of the EAD-UAV are continuous, indicating that the IBSA satisfies the constraints of the second-order continuity of the flight trajectory. It can be seen from Figure 19 that the attitude angles of the EAD-UAV obtained using IBSA are also continuous and within a reasonable range. As can be seen from Figure 20, the angle of attack and sideslip angle of the EAD-UAV are within the limited range. It can be seen from Figure 21 and Figure 22 that the flight trajectory obtained by IBSA optimization also satisfies the EAD thruster voltage and thrust constraints. The flight path determined by trajectory optimization using the IBSA is a second-order continuous curve; the flight time from the start to target 1 is 62.8919 s, that from target 1 to target 2 is 60.1779 s, and that from target 2 to target 3 is 62.7654 s; the calculation time of the optimization process is 6.3757 s.
Table A6 lists the flight times resulting from the mathematical simulation of the IBSA and BSA trajectory planning when Umax = [50:2:80] kV. The IBSA results are convergent within this range, which shows that the flight trajectory can be controlled by adjusting only the EAD-UAV voltage. With increasing Umax, the flight time decreases. A larger Umax makes the EAD-UAV have a higher flight speed and greater maneuverability. This result is in agreement with the expectations. The simulation using IBSA can converge to obtain the optimal trajectory with a Umax higher than 50 kV and a thrust limit higher than 5.2735 N. This indicates that the trajectory optimization of the EAD-UAV using IBSA can be realized with a thrust level similar to that in [9], which shows the feasibility of IBSA. As the BSA does not consider the subsequent optimization process when optimizing the trajectory of the current segment, it will leave a relatively poor initial value for the subsequent trajectory optimization problem, and the dynamic model of the EAD-UAV is relatively complex. Therefore, when using the BSA for trajectory optimization, only the first segment or first two segments of the three-segment flight trajectory converge, and the third segment flight trajectory diverges.

5. Conclusions

This paper presented the configuration of an EAD-UAV and the derivation of its attitude–path coupling dynamic equation. Based on this dynamic equation, the BSA method was used to optimize the three-dimensional flight trajectory between two points, and an IBSA algorithm was proposed to deal with the optimization of the multi-target flight trajectory rapidly when concerning the continuity of acceleration. The relationships among the boundary constraints, intermediate constraints, and Bezier basis function coefficients were deduced, and the continuous multi-target trajectory optimization problem was transformed into a single NLP problem that naturally satisfied the boundary condition and intermediate constraints. For the BSA used in the single-target scenario and IBSA, the simulation can converge with a Umax higher than 50 kV and a thrust limit higher than 5.2735 N. This level of thrust is similar to that in [9], and we can optimize the EAD thruster to reach this level, which indicates the feasibility of the BSA in single-target scenarios and IBSA for trajectory optimization of the EAD-UAV. The simulation showed that using the BSA to optimize the 3D trajectory of an EAD-UAV yielded results 1.14% different from the optimized performance index of GPM and a calculation time that was only 1.95% of that of GPM. Using the minimum energy consumption as the optimization goal, the average power was 1.4 kW, which is achievable. In the case of optimal time, the average power was four to eight times that in the case of the optimal energy consumption, leading to very high requirements for the battery. Therefore, balancing the energy consumption and flight time is very important when optimizing the flight trajectory of an EAD-UAV. Hence, the IBSA can overcome the poor convergence issue of the BSA under the continuous acceleration constraint for multi-target flight trajectories. For the EAD-UAV with the coupled dynamics, the IBSA can rapidly produce 3D trajectory optimization results.

Author Contributions

Conceptualization, M.H. and T.L.; methodology, M.H. and T.L.; software, M.H. and T.L.; validation, T.L., M.H. and N.Q.; formal analysis, M.H. and T.L.; investigation, M.H., T.L. and J.W.; resources, M.H.; data curation, T.L. and T.W.; writing—original draft preparation, M.H. and T.L.; writing—review and editing, T.L., Y.Z. and H.G.; visualization, M.H. and T.L.; supervision, N.Q. and M.H.; project administration, N.Q.; funding acquisition, N.Q. and M.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the National Natural Science Foundation of China (grant number 12272104) and in part by the National Natural Science Foundation of China (grant number U22B2013).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1. Airframe and EAD thruster parameters of the EAD-UAV.
Table A1. Airframe and EAD thruster parameters of the EAD-UAV.
Airframe parametersTotal mass (kg)2.6
Wingspan (m)5.14
Characteristic area (m2)4.8
Lift coefficient0.24
Drag coefficient0.03
Moment of inertia (kg·m2)Jx = 2.8
Jy = 0.4
Jz = 1.6
Jxy = 0.17
EAD thruster parametersRadius of emitting electrode(mm)0.1
Airfoil of collecting electrodeNACA0010
Gap between electrodes(mm)60
Span of electrode(m)3
Dimensionless constant C00.7
Ion mobility μ (m2·V−1·s−1)3 × 10−4
(cited from [51])
Number of electrode pairs in each thruster8
Thrust center distance (m)l1 = 0.1
l2 = 0.1
l3 = 0.2
Table A2. Parameters of the start and target in the single-target flight scenario.
Table A2. Parameters of the start and target in the single-target flight scenario.
Objective x (m) y (m) z (m) V (m/s) θ (°) ψ V (°) ϑ (°) ψ (°) γ (°) ω x (°/s) ω y (°/s) ω z (°/s)
Start0200500000000
Target1500220200500000000
Table A3. Parameters of the start and targets in the multi-target flight scenario.
Table A3. Parameters of the start and targets in the multi-target flight scenario.
Objective x (m) y (m) z (m) V (m/s) θ (°) ψ V (°) ϑ (°) ψ (°) γ (°) ω x (°/s) ω y (°/s) ω z (°/s)
Start0200500000000
Target 150012050
Target 21000120150
Target 31500220200500000000
Table A4. Flight and calculation times for single-target trajectory optimization using the BSA and GPM.
Table A4. Flight and calculation times for single-target trajectory optimization using the BSA and GPM.
Umax (kV)Fmax (N)BSAGPM
Ttotal (s)Calculation Time (s)Ttotal (s)Calculation Time (s)
505.2735231.21050.6726229.085032.8628
525.7436228.58990.7488226.343340.4004
546.2337225.95910.6525224.193732.7720
566.7437223.31790.6856221.097137.7848
587.2736220.67550.6389219.031631.0693
607.8234218.03380.7345215.689431.2783
628.3932215.39790.6979212.913331.8652
648.9830212.77280.6976210.334136.7595
669.5926210.16080.7502207.612738.3360
6810.2222207.56600.7256205.082334.9859
7010.8717204.99270.6866202.418041.3066
7211.5412202.42390.7748199.871337.7527
7412.2306199.89260.7923197.325540.7273
7612.9400197.41320.8059194.841838.5275
7813.6692194.94200.7460192.215745.8105
8014.4184192.50290.7824189.887449.3636
Table A5. Flight time, energy consumption, and average power of single-target trajectory optimization for JT and JEnergy.
Table A5. Flight time, energy consumption, and average power of single-target trajectory optimization for JT and JEnergy.
Umax (kV)Fmax (N)JTJEnergy
Ttotal (s)Energy Consumption
(W·h)
Average Power
(W)
Ttotal (s)Energy Consumption
(W·h)
Average Power
(W)
607.8234218.0338323.99825.3496 × 103267.8697106.66421.4335 × 103
628.3932215.3979345.61795.7764 × 103267.8685106.66371.4335 × 103
648.9830212.7728368.52846.2353 × 103267.8684106.66371.4335 × 103
669.5926210.1608391.75726.7107 × 103267.8683106.66371.4335 × 103
6810.2222207.5660416.71187.2274 × 103267.8686106.66381.4335 × 103
7010.8717204.9927442.79567.7762 × 103267.8681106.66361.4335 × 103
7211.5412202.4239470.23078.3628 × 103267.8712106.66481.4335 × 103
7412.2306199.8926498.19348.9723 × 103267.8683106.66371.4335 × 103
7612.9400197.4132527.24139.6147 × 103267.8700106.66431.4335 × 103
7813.6692194.9420557.53411.0296 × 104267.8678106.66351.4335 × 103
8014.4184192.5029589.32621.1021 × 104267.8711106.66481.4335 × 103
Table A6. Flight and calculation times for multi-target trajectory optimization using the IBSA and BSA.
Table A6. Flight and calculation times for multi-target trajectory optimization using the IBSA and BSA.
Umax (kV)IBSABSA
ΔT(1)
(s)
ΔT(2)
(s)
ΔT(3)
(s)
Ttotal
(s)
Calculation
Time (s)
ΔT(1)
(s)
ΔT(2)
(s)
ΔT(3)
(s)
Ttotal
(s)
Calculation
Time (s)
5076.413273.933976.5223226.86956.873575.5783invalidinvalidinvalidinvalid
5275.437772.903575.4856223.82685.829174.9793invalidinvalidinvalidinvalid
5474.617271.895274.5277221.04015.269574.0551invalidinvalidinvalidinvalid
5673.533870.966973.5909218.09176.743473.1152invalidinvalidinvalidinvalid
5872.812470.039472.6541215.50585.700772.1822invalidinvalidinvalidinvalid
6071.820069.087971.7320212.63995.955671.2523invalidinvalidinvalidinvalid
6270.771468.148770.7778209.69786.750570.3290invalidinvalidinvalidinvalid
6469.814567.222769.8679206.90517.041869.4099Invalid invalidinvalidinvalid
6668.929766.296968.9316204.15827.134368.4966invalidinvalidinvalidinvalid
6868.102565.402868.0002201.50566.627267.5913invalidinvalidinvalidinvalid
7067.294964.530567.1137198.93915.674866.6964invalidinvalidinvalidinvalid
7266.300763.626066.2121196.13876.793966.8113invalidinvalidinvalidinvalid
7465.345662.732265.3454193.42327.041264.9317invalidinvalidinvalidinvalid
7664.634561.889164.4766191.00026.100264.0640invalidinvalidinvalidinvalid
7863.771461.037963.6160188.42546.239263.2145invalidinvalidinvalidinvalid
8062.891960.177962.7654185.83526.375762.3733invalidinvalidinvalidinvalid

References

  1. Crowther, W.J.; Wilde, P.I.A.; Gill, K.; Michie, S.M. Towards integrated design of fluidic flight controls for a flapless aircraft. Aeronaut. J. 2009, 113, 699–713. [Google Scholar] [CrossRef]
  2. Wilde, P.I.A.; Buonanno, A.; Crowther, W.J.; Savvaris, A. Aircraft control using fluidic maneuver effectors. In Proceedings of the 26th AlAA Applied Aerodynamics Conference, Honolulu, HI, USA, 18–21 August 2008. [Google Scholar]
  3. Hoholis, G.; Steijl, R.; Badcock, K. Circulation control as a roll effector for unmanned combat aerial vehicles. J. Aircr. 2016, 53, 1875–1889. [Google Scholar] [CrossRef]
  4. Uluskan, S. Noncausal trajectory optimization for real-time range-only targetlocalization by multiple UAVs. Aerosp. Sci. Technol. 2020, 99, 105558. [Google Scholar] [CrossRef]
  5. Radmanesh, R.; Kumar, M.; French, D.; Casbeer, D. Towards a PDE-based large-scaledecentralized solution for path planning of UAVs in shared airspace. Aerosp. Sci. Technol. 2020, 105, 105965. [Google Scholar] [CrossRef]
  6. Yao, W.; Qi, N.; Wan, N.; Liu, Y. An iterative strategy for task assignment and path planning of distributed multiple unmanned aerial vehicles. Aerosp. Sci. Technol. 2019, 86, 455–464. [Google Scholar] [CrossRef]
  7. Liu, Y.; Wang, H.; Fan, J.; Wu, J.; Wu, T. Control-oriented UAV highly feasible trajectory planning: A deeplearning method. Aerosp. Sci. Technol. 2021, 110, 106435. [Google Scholar] [CrossRef]
  8. Sun, X.; Zhang, B.H.; Chai, R.Q.; Tsourdos, A.; Chai, S. UAV trajectory optimization using chance-constrained second-order cone programming. Aerosp. Sci. Technol. 2022, 121, 107283. [Google Scholar] [CrossRef]
  9. Xu, H.F.; He, Y.O.; Strobel, K.L. Flight of an aeroplane with solid-state propulsion. Nature 2018, 563, 532–539. [Google Scholar] [CrossRef]
  10. Drew, D.S.; Lambert, N.O.; Schindler, C.B. Toward controlled flight of the ionocraft: A flying microrobot using electrohydrodynamic thrust with onboard sensing and no moving parts. IEEE Robot. Autom. Lett. 2018, 3, 2807–2813. [Google Scholar] [CrossRef]
  11. Chattock, A.P.; Walker, W.E.; Dixon, E.H. On the specific velocities of ions in the discharge from points. Phil. Mag. Lett. 1901, 1, 79–98. [Google Scholar] [CrossRef]
  12. Stuetzer, O.M. Ion-drag pumps. J. Appl. Phys. 1960, 31, 136–146. [Google Scholar] [CrossRef]
  13. Brown, T.T. A Method of and an Apparatus or Machine for Producing Force or Motion. Great Britain Patent GB300311(A), 15 November 1928. [Google Scholar]
  14. Monrolin, N.; Ploouraboue, F.; Praud, O. Electrohydrodynamic thrust for in-atmosphere propulsion. AIAA J. 2017, 55, 4296–4305. [Google Scholar] [CrossRef]
  15. Gilmore, C.K.; Barrett, S.R. Electrohydrodynamic thrust density using positive corona-induced ionic winds for in-atmosphere propulsion. Proc. R. Soc. A Math. Phys. Eng. Sci. 2015, 471, 20140912. [Google Scholar] [CrossRef]
  16. Colas, D.F.; Ferret, A.; Pai, D.Z.; Lacoste, D.A.; Laux, C.O. Ionic wind generation by a wire-cylinder-plate corona discharge in air at atmospheric pressure. J. Appl. Phys. 2010, 108, 103306. [Google Scholar] [CrossRef]
  17. Kim, C.; Park, D.; Noh, K.C.; Hwang, J. Velocity and energy conversion efficiency characteristics of ionic wind generator in a multistage configuration. J. Electrost. 2010, 68, 36–41. [Google Scholar] [CrossRef]
  18. Miller, W.M. Force characterization of asymmetrical capacitor thrusters in air. AIAA J. 2009, 227, 293–327. [Google Scholar]
  19. Masuyama, K.; Barrett, S.R. On the Performance of Electrohydrodynamic Propulsion. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 2013, 469, 20120623. [Google Scholar] [CrossRef]
  20. Rickard, M.; Dunn-Rankin, D.; Weinberg, F.; Carleton, F. Characterization of ionic wind velocit. J. Electrost. 2005, 63, 711–716. [Google Scholar] [CrossRef]
  21. Rickard, M.; Dunn-Rankin, D.; Weinberg, F.; Carleton, F. Maximizing ion-driven gas flows. J. Electrost. 2006, 64, 368–376. [Google Scholar] [CrossRef]
  22. Wilson, J.; Perkins, H.D.; Thompson, W.K. An Investigation of Ionic Wind Propulsion; No. NASA/TM 2009-215822; NASA: Washington, DC, USA, 2009. [Google Scholar]
  23. Wynsberghe, E.; Turak, A. Station-keeping of a high-altitude balloon with electric propulsion and wireless power transmission: A concept study. Acta Astronaut. 2016, 128, 616–627. [Google Scholar] [CrossRef]
  24. Drew, D.S.; Pister, K.S. First take off of a flying micro-robot with no moving parts. In Proceedings of the 2017 International Conference on Manipulation, Automation and Robotics at Small Scales (MARSS), Piscataway, NJ, USA, 17–21 July 2017. [Google Scholar]
  25. He, Z.Z.; Li, P.F.; Wang, W.; Chen, X. Design of indoor unmanned airship propelled by ionic wind. J. Phys. Conf. Ser. 2021, 1748, 062011. [Google Scholar] [CrossRef]
  26. Coutinho, W.P.; Battarra, M.; Fliege, J. The unmanned aerial vehicle routing and trajectory optimisation problem, a taxonomic review. Comput. Ind. Eng. 2018, 120, 116–128. [Google Scholar] [CrossRef]
  27. Bollino, K.P.; Lewis, L.R.; Sekhavat, P.; Ross, I.M. Pseudospectral optimal control: A clear road for autonomous intelligent path planning. In Proceedings of the Infotech@ Aerospace 2007 Conference and Exhibit, Rohnert Park, CA, USA, 7–10 May 2007. [Google Scholar]
  28. Gong, Q.; Lewis, L.R.; Ross, I.M. Pseudospectral motion planning for autonomous vehicles. J. Guid. Control Dyn. 2009, 32, 1039–1045. [Google Scholar] [CrossRef]
  29. Bollino, K.P.; Lewis, L.R. Collision-free multi-UAV optimal path planning and cooperative control for tactical applications. In Proceedings of the Guidance, Navigation and Control Conference and Exhibit, Honolulu, HI, USA, 18–21 August 2008. [Google Scholar]
  30. Rogowski, K.; Maronski, R. Optimization of glider’s trajectory for given thermal conditions. Arch. Mech. Eng. 2011, 58, 11–25. [Google Scholar] [CrossRef]
  31. Pepy, R.; Herisse, B. An indirect method for optimal guidance of a glider. IFAC Proc. Vol. 2014, 47, 5097–5102. [Google Scholar] [CrossRef]
  32. Babel, L. Trajectory planning for unmanned aerial vehicles: A network optimization approach. Math. Methods Oper. Res. 2011, 74, 343–360. [Google Scholar] [CrossRef]
  33. Shanmugavel, M.; Tsourdos, A.; Zbikowski, R.; White, B.A.; Lechevin, N. A solution to simultaneous arrival of multiple UAVs using Pythagorean Hodograph curves. In Proceedings of the 2006 American Control Conference, Minneapolis, MN, USA, 14–16 June 2006. [Google Scholar]
  34. Subchan, S.; White, B.A.; Tsourdos, A.; Shanmugavel, M.; Zbikowski, R. Pythagorean Hodograph (PH) path planning for tracking airborne contaminant using sensor swarm. In Proceedings of the 2008 IEEE Instrumentation and Measurement Technology Conference, Minneapolis, MN, USA, 12–15 May 2008. [Google Scholar]
  35. Shah, M.A.; Aouf, N. 3D cooperative Pythagorean hodograph path planning and obstacle avoidance for multiple UAVs. In Proceedings of the 2010 IEEE 9th International Conference on Cyberntic Intelligent Systems, Reading, UK, 1–2 September 2010. [Google Scholar]
  36. Choe, R.; Puig-Navarro, J.; Cichella, V.; Xargay, E.; Hovakimyan, N. Cooperative trajectory generation using pythagorean hodograph bezier curves. J. Guid. Control Dyn. 2016, 39, 1744–1763. [Google Scholar] [CrossRef]
  37. Petropoulos, A.E.; Longuski, J.M. Shape-based algorithm for the automated design of low-thrust gravity assist trajectories. J. Spacecr. Rocket. 2004, 41, 787–796. [Google Scholar] [CrossRef]
  38. Jolly, K.G.; Sreerama, R.; Vijayakumar, R. A Bezier curve based path planning in a multi-agent robot soccer system without violating the acceleration limits. Robot. Auton. Syst. 2009, 57, 23–33. [Google Scholar] [CrossRef]
  39. Machmudah, A.; Parman, S.; Zainuddin, A. UAV Bezier curve maneuver planning using genetic algorithm. In Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation, Portland, OR, USA, 7–11 July 2010. [Google Scholar]
  40. Habib, Z.; Sakai, M. G2cubic transition between two circles with shape control. J. Comput. Appl. Math. 2009, 223, 134–144. [Google Scholar] [CrossRef]
  41. Connors, J.; Elkaim, G. Analysis of a spline based, obstacle avoiding path planning algorithm. In Proceedings of the VTC2007-Spring, IEEE 65th Vehicular Technology Conference, Dublin, Ireland, 22–25 April 2007. [Google Scholar]
  42. Huo, M.Y.; Yu, Z.; Liu, H.; Zhao, C.; Lin, T.; Song, Z.; Qi, N.M. Initial three-dimensional trajectory design for solar sails using Bezier shaping approach. IEEE Access 2019, 7, 150842–150850. [Google Scholar] [CrossRef]
  43. Huo, M.Y.; Mengali, G.; Quarta, A.; Qi, N.M. Electric sail trajectory design with Bezier curve-based shaping approach. Aerosp. Sci. Technol. 2019, 88, 126–135. [Google Scholar] [CrossRef]
  44. Chen, C.; He, Y.Q.; Bu, C.G.; Han, J.D.; Zhang, X. Quartic Bezier curve based trajectory generation for autonomous vehicles with curvature and velocity constraints. In Proceedings of the 2014 IEEE International Conference on Robotics and Automation, Hong Kong, China, 31 May–7 June 2014. [Google Scholar]
  45. Xiao, W.; Peng, J.; De, L.; Tao, S. Curvature continuous and bounded path planning for fixed-wing UAVs. Sensors 2017, 17, 2155. [Google Scholar]
  46. Yu, W.; Shuo, W.; Rui, W.; Min, T. Generation of temporal–spatial Bezier curve for simultaneous arrival of multiple unmanned vehicles. Inform. Sci. 2017, 418–419, 34–35. [Google Scholar] [CrossRef]
  47. Roth, J.R. Industrial Plasma Engineering Volume 1: Principles; Institute of Physics Publishing: Bristol, UK; Philadelphia, PA, USA, 1995; p. 286. [Google Scholar]
  48. Peek, F.W. Dielectric Phenomena in High Voltage Engineering; McGraw-Hill: New York, NY, USA, 1920; p. 57. [Google Scholar]
  49. Yu, Z.; Qi, N.M.; Huo, M.Y.; Fan, Z.C.; Yao, W.R. Fast Cooperative Trajectory Generation of Unmanned Aerial Vehicles Using a Bezier Curve-Based Shaping Method. IEEE Access 2022, 10, 1626–1636. [Google Scholar] [CrossRef]
  50. Yang, S.B.; Cui, T.; Hao, X.Y.; Yu, D.R. Trajectory optimization for a ramjet-powered vehicle in ascent phase via the Gauss pseudospectral method. Aerosp. Sci. Technol. 2017, 67, 88–95. [Google Scholar] [CrossRef]
  51. Moreau, E.; Benard, N.; Lan-Sun-Luk, J.D.; Chabriat, J.P. Electrohydrodynamic force produced by a wire-to-cylinder dc corona discharge in air at atmospheric pressure. J. Phys. D Appl. Phys. 2013, 46, 475204. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of the EAD propulsion principle.
Figure 1. Schematic diagram of the EAD propulsion principle.
Aerospace 10 00950 g001
Figure 2. EAD-UAV configuration with an orthogonal arrangement of multiple EAD thrusters.
Figure 2. EAD-UAV configuration with an orthogonal arrangement of multiple EAD thrusters.
Aerospace 10 00950 g002
Figure 3. Relationship between O x g y g z g and O x t y t z t .
Figure 3. Relationship between O x g y g z g and O x t y t z t .
Aerospace 10 00950 g003
Figure 4. Relationship between O x g y g z g and O x h y h z h .
Figure 4. Relationship between O x g y g z g and O x h y h z h .
Aerospace 10 00950 g004
Figure 5. Relationship between O x V y V z V and O x t y t z t .
Figure 5. Relationship between O x V y V z V and O x t y t z t .
Aerospace 10 00950 g005
Figure 6. Optimal single-target trajectory of an EAD-UAV for J = JT with Umax = 80 kV.
Figure 6. Optimal single-target trajectory of an EAD-UAV for J = JT with Umax = 80 kV.
Aerospace 10 00950 g006
Figure 7. Attitude of the EAD-UAV in the single-target optimal trajectory for J = JT with Umax = 80 kV.
Figure 7. Attitude of the EAD-UAV in the single-target optimal trajectory for J = JT with Umax = 80 kV.
Aerospace 10 00950 g007
Figure 8. Angle of attack and sideslip angle of EAD-UAV in the single-target optimal trajectory for J = JT with Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Figure 8. Angle of attack and sideslip angle of EAD-UAV in the single-target optimal trajectory for J = JT with Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Aerospace 10 00950 g008
Figure 9. Voltages of EAD-UAV thrusters with single-target optimal trajectories for J = JT obtained using Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Figure 9. Voltages of EAD-UAV thrusters with single-target optimal trajectories for J = JT obtained using Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Aerospace 10 00950 g009aAerospace 10 00950 g009b
Figure 10. Thrusts of EAD-UAV thrusters with single-target optimal trajectories for J = JT obtained using Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Figure 10. Thrusts of EAD-UAV thrusters with single-target optimal trajectories for J = JT obtained using Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Aerospace 10 00950 g010aAerospace 10 00950 g010b
Figure 11. Optimal single-target trajectory of an EAD-UAV for J = JEnergy with Umax = 60 kV and Umax = 80 kV.
Figure 11. Optimal single-target trajectory of an EAD-UAV for J = JEnergy with Umax = 60 kV and Umax = 80 kV.
Aerospace 10 00950 g011
Figure 12. Attitude of the EAD-UAV in the single-target optimal trajectory for J = JEnergy with Umax = 60 kV and Umax = 80 kV.
Figure 12. Attitude of the EAD-UAV in the single-target optimal trajectory for J = JEnergy with Umax = 60 kV and Umax = 80 kV.
Aerospace 10 00950 g012
Figure 13. Angle of attack and sideslip angle of EAD-UAV in the single-target optimal trajectory for J = JEnergy with Umax = 60 kV and Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Figure 13. Angle of attack and sideslip angle of EAD-UAV in the single-target optimal trajectory for J = JEnergy with Umax = 60 kV and Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Aerospace 10 00950 g013
Figure 14. Voltages of EAD-UAV thrusters with single-target optimal trajectories for J = JEnergy obtained using Umax = 60 kV and Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Figure 14. Voltages of EAD-UAV thrusters with single-target optimal trajectories for J = JEnergy obtained using Umax = 60 kV and Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Aerospace 10 00950 g014aAerospace 10 00950 g014b
Figure 15. Thrusts of EAD-UAV thrusters with single-target optimal trajectories for J = JEnergy obtained using Umax = 60 kV and Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Figure 15. Thrusts of EAD-UAV thrusters with single-target optimal trajectories for J = JEnergy obtained using Umax = 60 kV and Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Aerospace 10 00950 g015aAerospace 10 00950 g015b
Figure 16. Thrusts of EAD-UAV thrusters with multi-target optimal trajectory obtained using Umax = 80 kV.
Figure 16. Thrusts of EAD-UAV thrusters with multi-target optimal trajectory obtained using Umax = 80 kV.
Aerospace 10 00950 g016
Figure 17. Velocity of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Figure 17. Velocity of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Aerospace 10 00950 g017
Figure 18. Acceleration of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Figure 18. Acceleration of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Aerospace 10 00950 g018
Figure 19. Attitude of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Figure 19. Attitude of an EAD-UAV with the multi-target optimal trajectory obtained using Umax = 80 kV.
Aerospace 10 00950 g019
Figure 20. Angle of attack and sideslip angle of EAD-UAV in the multi-target optimal trajectory with Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Figure 20. Angle of attack and sideslip angle of EAD-UAV in the multi-target optimal trajectory with Umax = 80 kV. (a) Angle of attack, (b) sideslip angle.
Aerospace 10 00950 g020
Figure 21. Voltages of EAD-UAV thrusters with multi-target optimal trajectory obtained using Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Figure 21. Voltages of EAD-UAV thrusters with multi-target optimal trajectory obtained using Umax = 80 kV. (a) Voltage of thruster 1, (b) voltage of thruster 2, (c) voltage of thruster 3, (d) voltage of thruster 4, (e) voltage of thruster 5, (f) voltage of thruster 6.
Aerospace 10 00950 g021
Figure 22. Thrusts of EAD-UAV thrusters for multi-target optimal trajectory with Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Figure 22. Thrusts of EAD-UAV thrusters for multi-target optimal trajectory with Umax = 80 kV. (a) Thrust of thruster 1, (b) thrust of thruster 2, (c) thrust of thruster 3, (d) thrust of thruster 4, (e) thrust of thruster 5, (f) thrust of thruster 6.
Aerospace 10 00950 g022aAerospace 10 00950 g022b
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lin, T.; Huo, M.; Qi, N.; Wang, J.; Wang, T.; Gu, H.; Zhang, Y. Coupling Dynamics and Three-Dimensional Trajectory Optimization of an Unmanned Aerial Vehicle Propelled by Electroaerodynamic Thrusters. Aerospace 2023, 10, 950. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace10110950

AMA Style

Lin T, Huo M, Qi N, Wang J, Wang T, Gu H, Zhang Y. Coupling Dynamics and Three-Dimensional Trajectory Optimization of an Unmanned Aerial Vehicle Propelled by Electroaerodynamic Thrusters. Aerospace. 2023; 10(11):950. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace10110950

Chicago/Turabian Style

Lin, Tong, Mingying Huo, Naiming Qi, Jianfeng Wang, Tianchen Wang, Haopeng Gu, and Yiming Zhang. 2023. "Coupling Dynamics and Three-Dimensional Trajectory Optimization of an Unmanned Aerial Vehicle Propelled by Electroaerodynamic Thrusters" Aerospace 10, no. 11: 950. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace10110950

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop