Next Article in Journal
Multi-Layer Fault-Tolerant Robust Filter for Integrated Navigation in Launch Inertial Coordinate System
Next Article in Special Issue
Fault-Tolerant Control for Hexacopter UAV Using Adaptive Algorithm with Severe Faults
Previous Article in Journal
Space Environment Evaluation and Low-Earth-Orbit Demonstration of a Communication Component with a Commercial Transceiver Integrated Circuit
Previous Article in Special Issue
The Flying and Adhesion Robot Based on Approach and Vacuum
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Aggressive Constrained Trajectory Synthesis and Control for Multi-Copters

by
Tsung-Liang Liu
and
Kamesh Subbarao
*,†
Department of Mechanical and Aerospace Engineering, The University of Texas at Arlington, 500 W. First St, Box 19018, 211 Woolf Hall, Arlington, TX 76019, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 11 April 2022 / Revised: 11 May 2022 / Accepted: 17 May 2022 / Published: 24 May 2022
(This article belongs to the Special Issue Applications of Drones)

Abstract

:
In this paper, we propose a novel time and control effort optimal aggressive trajectory synthesis and control design methodology. The trajectory synthesis is a modified minimum snap design but with specific position and orientation constraints on a multi-copter, such as flying through tight spaces (windows) at specific orientations. The paper also introduces a means to stitch together multiple flight segments, enforce smoothness, and minimize segment times as well as the overall time, thereby resulting in very aggressive and feasible trajectories. A novel analysis for a specific scenario when no yaw angle specifications are provided is conducted, wherein a trade-off results in additional aggressiveness. The control algorithms to follow these trajectories are based on an inverse dynamics approach. Several candidate high-fidelity simulations are performed to verify the effectiveness of the proposed approach.

1. Introduction

Multi-copters, such as quadcopters, hexacopters, octocopters, and others have been very popular in recent years mainly because of their availability, mobility, agility, and flexibility. They are a great platform for control experiments and various applications. Their characteristics also make them an attractive choice for high-speed aerial navigation through complex environments. Various researches related to aggressive trajectory generation and aggressive maneuver tracking have been conducted. Some recent examples include [1], in which the authors proposed a novel control law for accurate tracking of aggressive quadcopter trajectories. In [2] the authors presented a framework to do optimal time allocation for quadcopter trajectory generation. In [3] the authors addressed the problem of performing aggressive quadcopter maneuvers that are attitude-constrained.
Our main contribution in this paper is the development of a synthesis framework that fully utilizes the available dynamic capability of the multi-copter when performing the most aggressive maneuver. We pursue aggressiveness because time is a critical issue in the given scenario. We adopt aggressive maneuvers because the environment is complicated and sometimes very specific maneuvers are needed in order to satisfy path and vehicle state constraints, given the environment model. However, the level of aggressiveness one can achieve depends on the optimality of the trajectory plan and the dynamic capability of the vehicle. It will be valuable if one can generate a trajectory that fully utilizes the dynamic capability of the vehicle while accommodating the maneuvers required at specific locations. Therefore, in this research, we assume that the dynamic capability of the multi-copter is limited by its maximum rotor thrust and define the aggressiveness accordingly.
For multi-copter trajectory generation and motion control, ideas such as multi-segment polynomial, minimum derivative optimization, and differential flatness inverse dynamics analysis are widely adopted [4,5,6]. To achieve an aggressive trajectory, methods such as segment time allocation [7] and spatial-temporal trajectory [2] have been proposed. However, we do not yet have a synthetic architecture to find the most aggressive trajectory based on the waypoints, the requirement for vehicle heading and maneuvering at waypoints, and the constraints of rotor thrust for a multi-coper. We would further like to emphasize that the contribution of this research lies in the development of a complete synthesis to fulfill this gap. The constraint on vehicle maneuvering at waypoints is posed as a window passing through problem and solved through the trajectory optimization as trajectory derivative constraints. The maximum rotor thrust for the trajectory is obtained through the inverse dynamics analysis. As the optimized aggressive trajectory is obtained, a geometric controller is used to perform the polynomial trajectory tracking and verify the feasibility of the trajectory and the accuracy of the inverse dynamics analysis. Further, a novel method for yaw trajectory optimization is also proposed to further improve the aggressive performance in the case when there are no specific requirements on the vehicle heading.
The rest of the paper is organized as follows. Section 2 describes the problem and the setup. Section 3 summarizes the mathematical modeling providing the basic governing equations. Section 4 discusses the solution methodology in detail, followed by the numerical results in Section 5. Section 6 summarizes the performance of the geometric controller with Section 7 providing details on the modification to the trajectory should there be no heading constraints. Finally, Section 8 provides concluding remarks and summarizes future work.

2. Problem Description and Setup

While aggressive flight using multi-copters is widely discussed/mentioned, there is not much effort put toward defining aggressiveness and finding the most aggressive trajectory with complex maneuver requirements for a given vehicle model. Therefore, in this paper, the problem is addressed as: Given waypoint coordinates, heading angles, and some vehicle velocity and attitude constraints at these waypoints, how can we find a feasible trajectory that fully utilizes the dynamic capability of a multi-copter to achieve the optimal aggressiveness? A formal definition of aggressiveness will be provided in the latter sections.

3. Mathematical Model Description

A generic model for quadcopters is used in this paper. Note that the proposed approach can be applied to different configurations of multi-copters by modifying the vehicle force and moment model according to the desired rotor configuration and control allocation method. The coordinate systems and forces and moments generated by the rotors are shown in Figure 1.
The body frame, B, is attached to the center of mass of the quadcopter and rotor 1 is on the positive X B -axis. The gravitational acceleration g is in the Z I direction of the inertial frame, I. Euler angles roll ϕ , pitch θ , and yaw ψ are used to define orientation from inertial frame to body frame. Note, Z X Y rotation order is used here, and therefore the rotation matrix for transforming coordinates from B to I is given by:
R I B = cos ψ cos θ sin ϕ sin ψ sin θ cos ϕ sin ψ cos ψ sin θ + cos θ sin ϕ sin ψ cos θ sin ψ + cos ψ sin ϕ sin θ cos ϕ cos ψ sin ψ sin θ cos ψ cos θ sin ϕ cos ϕ sin θ sin ϕ cos ϕ cos θ
The position vector of the quadcopter in the inertial frame is denoted by r . With the gravity force acting in the Z I direction and the forces of the rotors acting in the Z B direction, the equation governing the acceleration of the quadcopter with respect to inertial frame is given by:
m r ¨ = 0 0 m g + R I B 0 0 F 1 + F 2 + F 3 + F 4
With p, q, and r denoting the components of angular velocity of the quadcopter in the body frame, the rotational kinematics equation is given by:
p q r = cos θ 0 cos ϕ sin θ 0 1 sin ϕ sin θ 0 cos ϕ cos θ ϕ ˙ θ ˙ ψ ˙
Assuming rotors 1 and 3 rotate in the Z B direction while 2 and 4 rotate in the Z B direction, M 1 and M 3 act in the Z B direction while M 2 and M 4 act in the Z B direction since the moment produced by the rotor is opposite to the direction of rotation of the blade. With I denoting the moment of inertia of the quadcopter referenced to the center of mass and L denoting the distance from the axis of rotation of the rotors to the center of the quadcopter, the rotational dynamics equation is given by:
I p ˙ q ˙ r ˙ = L ( F 2 F 4 ) L ( F 3 F 1 ) M 1 M 2 + M 3 M 4 p q r × I p q r
Define the input u = [ u 1 u 2 u 3 u 4 ] T wherein u 1 is the total force from the rotors and u 2 , u 3 and u 4 are the moments about X B , Y B and Z B axes. Following [6], we assume that the force and moment produced by the ith rotor are proportional to the square of its rotational speed as:
F i = k F ω i 2 , M i = k M ω i 2
The relationship between the input and the angular speed of the rotors can be represented as:
u 1 u 2 u 3 u 4 = F 1 + F 2 + F 3 + F 4 L ( F 2 F 4 ) L ( F 3 F 1 ) M 1 M 2 + M 3 M 4 = k F k F k F k F 0 k F L 0 k F L k F L 0 k F L 0 k M k M k M k M ω 1 2 ω 2 2 ω 3 2 ω 4 2

4. Solution Methodology

We use multi-segment polynomials to generate the trajectory and assure its smoothness by having trajectory derivative continuities. To optimize the polynomial coefficients, we form the quadratic problem with the cost on the trajectory derivative which is related to the actuator input of the vehicle. We accommodate waypoint velocity and attitude requirements in the optimization problem as equality constraints. To optimize the segment time between waypoints, we introduce an augmented cost on both trajectory derivative and total time. We connect the polynomial trajectory and actual quadcopter dynamics and rotor inputs by performing inverse dynamics analysis. We finally define aggressiveness and find the optimal trajectory with desired aggressiveness by tracking the desired maximum rotor force needed for the optimized trajectory. The whole methodology and the solution procedure are illustrated below in Figure 2.

4.1. Multi-Segment Polynomial Trajectory Optimization

Multi-segment polynomials are used to generate the trajectory. For each segment connecting one and the subsequent waypoint, independent polynomials are used to represent the quadcopter states x, y, z, and ψ (yaw angle). Each polynomial segment is represented as:
P ( t ) = p N t N + p N 1 t N 1 + + p 0 = i = 0 N p i t i
The cost function on the integral of the quadratic of the trajectory derivatives is:
J = 0 T [ c 0 P ( t ) 2 + c 1 P ( t ) 2 + c 2 P ( t ) 2 + + c N P ( N ) ( t ) 2 ] d t
where T is the flight time for the trajectory segment and c r is user-specified penalty on the r t h derivative of the trajectory. This function can be written in matrix form as:
J = p ¯ T Q p ¯
where p ¯ = [ p 0 p 1 p N ] T and Q is the cost matrix. Following the formulation in [4], Q can be constructed as:
for r = 0 , 1 , N , i = 0 , 1 , N , l = 0 , 1 , N Q r ( i + 1 , l + 1 ) = [ m = 0 r 1 ( i m ) ( l m ) ] T i + l 2 r + 1 i + l 2 r + 1 if i r l r 0 if i < r l < r
Q = r = 0 N c r Q r
In this paper, to minimize the input needed to achieve optimal aggressiveness, we adopted the choice in [5] which is to minimize the snap (4th derivative) in trajectory x, y, and z while minimizing the acceleration (2nd derivative) in trajectory ψ . Therefore, we use c 4 = 1 while all other coefficients are set to zero for x, y, and z polynomial, and only c 2 = 1 for yaw polynomial. The constraints on the derivatives on the endpoints of a polynomial segment can be imposed as a linear function of the coefficients:
A p ¯ = b A = A 0 A T , b = b 0 b T
where A is constructed by evaluating the components in the derivative formulations of the polynomial at t = 0 and t = T corresponding to the appropriate coefficients as:
A 0 r n = m = 0 r 1 ( n m ) if r = n 0 if r n
A T r n = [ m = 0 r 1 ( n m ) ] T n r if r n 0 if r > n
We use the 0th order derivative constraint to specify the waypoint position. Higher order derivatives can be used to specify desired waypoint velocity, acceleration, etc., e.g., to enforce that the quadcopter starts from rest at the beginning of a trajectory. If not specified, these derivatives are subject to minimization of the cost function. Having assembled Q , A , and b , the quadratic problem can be written below. There are methods to solve such a standard equality constrained QP [8]. In this research, the Matlab solver “quadprog” is used.
min p ¯ p ¯ T Q p ¯ s . t . A p ¯ b = 0
For M polynomial segments, the joint optimization can be composed by concatenating their cost matrices in a block-diagonal fashion as:
J joint = p ¯ 1 p ¯ M T Q 1 ( T 1 ) 0 0 0 0 0 0 Q M ( T M ) p ¯ 1 p ¯ M
The derivative constraints can also be concatenated in a block-diagonal fashion as:
A 1 ( T 1 ) 0 0 0 0 0 0 A M ( T M ) p ¯ 1 p ¯ M = A der p ¯ 1 p ¯ M = b 1 b M
To have a feasible smooth trajectory, the derivatives should be continuous between segments. The continuity constraints must be imposed to ensure that the derivatives at the end of the ith segment match the derivatives at the beginning of the ( i + 1 ) th segment:
A T , i p ¯ i = A 0 , i + 1 p ¯ i + 1
A T , 1 A 0 , 2 0 0 0 0 0 A T , 2 A 0 , 3 0 0 0 0 0 0 0 A T , M 1 A 0 , M p ¯ 1 p ¯ 2 p ¯ 3 p ¯ M 1 p ¯ M = A con p ¯ 1 p ¯ 2 p ¯ 3 p ¯ M 1 p ¯ M = 0 0 0 0 0
These constraints can be compiled into a single set of linear equality constraints for the joint optimization problem:
A der A con p ¯ 1 p ¯ M = b 1 b M 0 0

4.2. Constraints on Vehicle Velocity and Attitude at Waypoints

In scenarios such that the quadcopter needs to pass through some narrow gap (e.g., a window) as shown in Figure 3, we can utilize the derivative constraints to achieve the objective.
Assume a waypoint at the center of the window and the window orientation is defined with Z X Y Euler angles roll ϕ , pitch θ , and yaw ψ . The window forward vector in inertial frame can be obtained as:
w F = R I B 1 0 0 = cos ψ cos θ sin ϕ sin ψ sin θ cos θ sin ψ + cos ψ sin ϕ sin θ cos ϕ sin θ = W F x W F y W F z
The window upward vector in inertial frame can be obtained as:
w U = R I B 0 0 1 = cos ψ sin θ + cos θ sin ϕ sin ψ sin ψ sin θ cos ψ cos θ sin ϕ cos ϕ cos θ = W U x W U y W U z
To fly through the window safely, the velocity vector at the waypoint should be aligned with the window forward vector. Therefore, the cross product of the vectors should be zero, and the constraints for the waypoint 1st order derivatives can be composed as:
v × w F = x ˙ y ˙ z ˙ × W F x W F y W F z = y ˙ W F z z ˙ W F y z ˙ W F x x ˙ W F z x ˙ W F y y ˙ W F x = 0 0 0
From the equation of acceleration, we can deduce that z B in inertial frame is in the direction of the vector [ x ¨ y ¨ z ¨ + g ] T :
m x ¨ y ¨ z ¨ = 0 0 m g + R I B 0 0 1 F m x ¨ y ¨ z ¨ + g = z B · F
For safely flying through the window, the z B vector of the quadcopter should be aligned with the window upward vector. Therefore, the constraints for the waypoint 2nd order derivatives can be composed as:
x ¨ y ¨ z ¨ + g × W U x W U y W U z = y ¨ W U z ( z ¨ + g ) W U y ( z ¨ + g ) W U x x ¨ W U z x ¨ W U y y ¨ W U x = 0 0 0
y ¨ W U z z ¨ W U y z ¨ W U x x ¨ W U z x ¨ W U y y ¨ W U x = g W U y g W U x 0
For these constraints, the relationship between x, y, and z derivatives needs to be specified. Therefore, the joint optimization of x, y, and z should be composed by concatenating their cost and constraint matrices in a block-diagonal fashion as:
J x y z = p ¯ ¯ x p ¯ ¯ y p ¯ ¯ z T Q x 0 0 0 Q y 0 0 0 Q z p ¯ ¯ x p ¯ ¯ y p ¯ ¯ z
A x 0 0 0 A y 0 0 0 A z p ¯ ¯ x p ¯ ¯ y p ¯ ¯ z = b x b y b z
where p ¯ ¯ x = [ p ¯ x 1 p ¯ x M ] T , p ¯ ¯ y = [ p ¯ y 1 p ¯ y M ] T and p ¯ ¯ z = [ p ¯ z 1 p ¯ z M ] T . With this joint optimization setup, the constraints for a window passage on waypoint s in an M-segment trajectory can be implemented as:
0 W F z v 0 , s W F y v 0 , s W F z v 0 , s 0 W F x v 0 , s W F y v 0 , s W F x v 0 , s 0 0 W U z a 0 , s W U y a 0 , s W U z a 0 , s 0 W U x a 0 , s W U y a 0 , s W U x a 0 , s 0 p ¯ ¯ x p ¯ ¯ y p ¯ ¯ z = 0 0 0 g W U y g W U x 0
v 0 , s = [ 0 × ( n + 1 ) ] × ( s 1 ) v 0 [ 0 × ( n + 1 ) ] × ( M s ) a 0 , s = [ 0 × ( n + 1 ) ] × ( s 1 ) a 0 [ 0 × ( n + 1 ) ] × ( M s )
For a segment polynomial P ( t ) = p 0 + p 1 t 1 + p 2 t 2 + p 3 t 3 + + p N t N , the 1st order derivative (velocity) at t = 0 is v 0 p ¯ and v 0 = [ 0 1 0 0 0 ] . The 2nd order derivative (acceleration) at t = 0 is a 0 p ¯ and a 0 = [ 0 0 2 0 0 ] .

4.3. Segment Time Optimization

The above optimization finds the optimal trajectory polynomials for the given segment times. However, in general cases, we do not specify segment times, and instead, we would like to also optimize them to achieve better aggressiveness. To optimize segment times, we modify the cost function in the form:
J T = J o r g + c T i = 1 M T i
where J o r g combines the original cost on x, y, z, and ψ while c T is a penalty on total time. The cost on x, y, z is the integral of square of derivatives on distance while the cost on ψ is on angle. To combine the costs, two coefficients μ r and μ ψ are introduced to non-dimensionalize the costs.
J o r g = μ r J x y z + μ ψ J ψ
In this research, we minimize the 4th derivative in trajectory x, y, and z while minimizing the 2nd derivative in trajectory yaw angle. Therefore, the above equation can be represented as:
J o r g = μ r 0 T f d 4 x d t 4 2 + d 4 y d t 4 2 + d 4 z d t 4 2 d t + μ ψ 0 T f d 2 ψ d t 2 2 d t
where μ r and μ ψ are defined as:
μ r = 1 max d 4 x d t 4 , d 4 y d t 4 , d 4 z d t 4 2
μ ψ = 1 max d 2 ψ d t 2 2
We find μ r and μ ψ in the first iteration with initial segment times and then use them as constants throughout the optimization process. With the cost function defined, we perturb each segment time by some δ t to obtain the gradient of the cost function with respect to each segment time. This is then used in a gradient descent method to find the time allocation for the minimum cost iteratively.
T ¯ = T 1 T M T
i J T = J T T 1 T i + δ t T M T J T ( T ¯ ) δ t
T i n e w = T i + α Δ T i , Δ T i = i J T i J T
Since this is a high-dimensional problem with a complex cost function, the step size α can easily become too small or too large during the iterations and lead to slow convergence or even divergence. For numerical efficiency, stability, and convergence, we use the backtracking line search method [9] to find a suitable step size α in every iteration.
while J ( T ¯ + α Δ T ¯ ) > J ( T ¯ ) + ϵ α J T T Δ T ¯ , α : = β α
With β = 0.5 and ϵ = 0.0001 , we have fast and stable convergence in our test cases. Figure 4 shows the result of segment time optimization for a simple 2D 4-waypoint scenario. The initial segment time is T ¯ = [ 1 3 1 ] T s, and we optimize it with two different c T values, 50 and 5000. The result trajectories look identical because they have similar segment time distribution while the total time turned out to be 6.28 and 3.53 seconds, respectively. It is observed that, with this segment time optimization process, the optimal segment time distribution will be found to minimize the cost on integral of the square of trajectory derivatives, while the optimal total time is found based on the value of c T . The larger c T used, the smaller the total time will be, which makes the trajectory more aggressive.

4.4. Inverse Dynamics

To connect the polynomial trajectory optimization framework and actual quadcopter dynamics, we perform the inverse dynamics analysis. The differential flatness method is widely adopted in inverse dynamics analysis for multi-copters [6]. In this process, we will find all the states and inputs of the quadcopter according to the trajectory x, y, z, ψ , and their derivatives. For orientation ϕ and θ , first from the equation of acceleration we have:
z B = x ¨ y ¨ z ¨ + g T x ¨ y ¨ z ¨ + g , u 1 = m x ¨ y ¨ z ¨ + g
Assume x C is the vector obtained by rotating x I around z I by yaw angle ψ :
x C = cos ψ sin ψ 0 T
We can determine x B and y B by:
y B = z B × x C z B × x C , x B = y B × z B
With vehicle body frame defined, we can determine the rotation matrix and roll and pitch angles by:
R I B = x B y B z B
ϕ = arcsin [ R I B ( 3 , 2 ) ] , θ = arctan R I B ( 3 , 1 ) R I B ( 3 , 3 )
For angular velocity p, q, and r, first we take the 1st derivative of the equation of acceleration:
m a ˙ = u ˙ 1 z B + ω B × u 1 z B
Substituting u ˙ 1 = z B · m a ˙ we have:
ω B × z B = m u 1 [ a ˙ ( z B · a ˙ ) z B ]
With a ˙ = [ x y z ] , the RHS is known. With x B , y B , and z B being unit vectors, ω B × z B can be considered as the projection of ω B onto the x B y B plane with 90 shift. Therefore, the body angular velocities p and q can be determined as:
p = ( ω B × z B ) · y B , q = ( ω B × z B ) · x B
From the rotational kinematics equation, we have:
ψ ˙ = sin θ cos ϕ p + cos θ cos ϕ r r = ψ ˙ + sin θ cos ϕ p cos ϕ cos θ
With p, q, and r solved, we have ω B = R I B [ p q r ] T , and θ ˙ and ϕ ˙ can also be obtained from the inversion of the rotational kinematics equation. For angular acceleration p ˙ , q ˙ and r ˙ , first we take the 2nd derivative of the equation of acceleration:
m a ¨ = u ¨ 1 z B + 2 ω B × u ˙ 1 z B + ω B × ω B × u 1 z B + ω ˙ B × u 1 z B ω ˙ B × z B = ( m a ¨ u ¨ 1 z B 2 ω B × u ˙ 1 z B ω B × ω B × u 1 z B ) / u 1
With a ¨ = [ x y z ] , u ¨ 1 = z B · ( m a ¨ ω B × ω B × u 1 z B ) and u ˙ 1 = z B · m a ˙ , the RHS is known and the body angular accelerations p ˙ and q ˙ can be determined as:
p ˙ = ( ω ˙ B × z B ) · y B , q ˙ = ( ω ˙ B × z B ) · x B
Taking derivative of previous ψ ˙ equation, we have:
ψ ¨ = sin θ cos ϕ p ˙ + cos θ cos ϕ r ˙ p d d t sin θ cos ϕ + r d d t cos θ cos ϕ r ˙ = cos ϕ cos θ ψ ¨ + sin θ cos ϕ p ˙ + p cos θ cos ϕ θ ˙ + sin θ sin ϕ ϕ ˙ cos 2 ϕ r cos θ sin ϕ ϕ ˙ sin θ cos ϕ θ ˙ cos 2 ϕ
Next we determine the moment inputs from obtained angular velocity and acceleration by the rotational dynamics equation:
u 2 u 3 u 4 = I p ˙ q ˙ r ˙ + p q r × I p q r
Now we have found all the states and inputs of the quadcopter derived from the trajectory x, y, z, ψ , and their derivatives. We can further determine the angular speed and force produced of each rotor by:
ω 1 2 ω 2 2 ω 3 2 ω 4 2 = k F k F k F k F 0 k F L 0 k F L k F L 0 k F L 0 k M k M k M k M 1 u 1 u 2 u 3 u 4 , F i = k F ω i 2

4.5. Max Force Tracking and Aggressiveness Defined

We find the maximum rotor force needed for the trajectory through the inverse dynamics analysis and observe how the maximum force and the trajectory total time vary with different c T values. Figure 5 shows the analysis result for a 3D 4-waypoint scenario. It is observed that, with the increment in c T , the total time decreases smoothly while the maximum force rises smoothly with increasing slope. With this relationship, we can track a particular total time or maximum force for a given scenario by the gradient descent method through c T .
With optimal aggressiveness, the quadcopter should finish the tasks in the shortest time possible by means of its dynamic capability. Generally speaking, the dynamic capability of a quadcopter is restricted by the maximum thrust of its rotors. Therefore, we define aggressiveness as the percentage of excess thrust required for the optimized trajectory.
Aggressiveness = F m a x _ r e q F h o v e r F m a x F h o v e r · 100 %
where F m a x _ r e q is the maximum rotor force required for the trajectory, F h o v e r is the rotor force for steady hover, and F m a x is the maximum thrust of the rotor. One might conserve the aggressiveness for safety reasons. Given waypoints x, y, z, ψ , and the desired aggressiveness, the whole methodology to find the optimal trajectory and segment time allocation for a given quadcopter model is summarized in Figure 2.

5. Numerical Results

The quadcopter parameters we use in this research are from [10] and tabulated below in Table 1.
We test the optimization framework with a four-waypoint scenario. The waypoint settings are tabulated below (yaw angles at waypoints 2 and 3 are not specified) in Table 2.
We specify the mission to be from rest to rest, therefore the velocity and acceleration at the first and the last waypoint are constrained to be zero. There are two narrow windows to pass through in this scenario. The window settings are tabulated below in Table 3.
We specify the aggressiveness to be 80%. For this quadcopter model the rotor force for steady hover is F h o v e r = 2.5 N. Assuming F m a x = 3.75 N, we have the maximum rotor force we can use for the trajectory as F m a x _ r e q = 3.5 N. With initial guess of T ¯ = [ 5 5 5 ] T s and c T = 100 , we use 10th order segment polynomials and have derivative continuity constraints on up to 6th order derivative (Pop). Figure 6 shows the scenario and the optimal trajectory.
Figure 7, Figure 8 and Figure 9 show the result of the inverse dynamics analysis from the trajectory with respect to the particular quadcopter model. The magenta dashed lines mark the times of waypoint passages. It can be noticed that the quadcopter had a spike in angular acceleration and angular rate around waypoint 3 because it was making the maneuver to pass through the highly tilted window 2.
Figure 10 and Figure 11 show the control input and the rotor thrust along the trajectory. From the plot, we confirmed that the maximum rotor thrust used is 3.5 N as expected.
Figure 12 shows the quadcopter on the trajectory with a fixed time interval of 0.33 s. It can be observed that the quadcopter had the highest speed in segment 2 so it can maneuver through the highly tilted window 2. The successful passage through the narrow windows can be visually confirmed. If we further check the data at waypoint 3, we have the window forward and upward vectors as w F = [ 0.939 0.342 0 ] T and w F = [ 0.171 0.469 0.866 ] T . Additionally, the quadcopter velocity and thrust vectors (aligned with z B ) are v = [ 1.777 0.6470 0 ] T and f = [ 1.454 3.995 7.364 ] T . The velocity and thrust vectors are confirmed to be aligned with the window forward and upward vectors. Finally, the optimal time allocation was found as T ¯ o p t = [ 1.29 2.29 2.62 ] T s.

6. Geometric Control and Simulation Result

To perform the polynomial trajectory tracking and verify the result of the aggressive trajectory optimization, we adopted the geometric controller proposed by Lee et al. in [11] and implemented it in MATLAB/Simulink. The controller is constructed in two parts, trajectory tracking and attitude tracking. The control input of total force f is obtained by the trajectory tracking part with:
f = ( k x e x k v e v m g e 3 + m x ¨ d ) · R e 3
The control input of the moments is obtained by the attitude tracking part with:
M = k R e R k ω e ω + ω × I ω I ( ω ^ R T R d ω d R T R d ω ˙ d )
where vector e 3 defines the direction of gravity and R is the rotation matrix from body to inertial frame. The tracking errors for position e x , velocity e v , attitude e R , and angular velocity e ω are defined as:
e x = x x d e v = v v d
e R = 1 2 ( R d T R R T R d ) e ω = ω R T R d ω d
where the h a t m a p ^ : R 3 SO ( 3 ) is defined by the condition that x ^ y = x × y for all x , y R 3 . Additionally, the v e e m a p : SO ( 3 ) R 3 is the inverse of the hat map. At a given moment, x , v , ω , and R represent the position, linear velocity, body angular velocity, and rotation matrix of the vehicle. The corresponding desired position x d and velocity v d are captured from the polynomial trajectory. The desired body z axis can be obtained as:
z d = k x e x k v e v m g e 3 + m x ¨ d k x e x k v e v m g e 3 + m x ¨ d
With the desired body z axis, yaw angle, and derivatives of the polynomial trajectory available, we can find the desired rotation matrix R d , angular velocity ω d , and angular acceleration ω ˙ d by following the process discussed previously in inverse dynamics. The simulation is implemented with the following equations of motion:
x ˙ = v m v ˙ = m g e 3 f R e 3
R ˙ = R ω ^ I ω ˙ + ω × I ω = M
With 10 % error in estimated inertia, ± 1 N noise added to f and ± 0.003 Nm noise added to M (around 10 % of the maximum input used) as control input disturbance, we have the simulation result shown in figures below.
From the simulation result we can see that, despite the existence of inertia estimation error and control input disturbance, the geometric control had a good performance in tracking the polynomial trajectory. Position and yaw angle are well tracked and roll and pitch angles turned out to be very close to the inverse dynamics estimation as shown in Figure 13 and Figure 14. Therefore, the successful passage through the narrow windows is confirmed. Figure 15 and Figure 16 show that the linear acceleration is well tracked and the angular acceleration is very close to the estimation. From Figure 17, the control input used to track this polynomial trajectory is verified to be very close to the inverse dynamics estimation. The maximum rotor thrust used in the simulation is also confirmed to be around the maximum rotor force (3.5 N) requirement we specified in the trajectory generation as shown in Figure 18.

7. Yaw Trajectory Optimization

If we do not have a specific requirement on yaw angle, can we optimize the yaw trajectory to achieve better aggressive performance? In other words, by finding the yaw trajectory such that no M z control moment input is required to track the aggressive polynomial trajectory, we can further reduce the maximum rotor force needed and thereby track the aggressive trajectory with less motor force or achieve a faster trajectory with the same aggressiveness specified. From the equations of motion:
M x M y M z = I p ˙ q ˙ r ˙ + p q r × I p q r
Because of the symmetry, the moment of inertia of the quadcopter is assumed as:
I = I x x 0 0 0 I y y 0 0 0 I z z , and I x x = I y y
M z = I z z r ˙ + p q ( I y y I x x ) = I z z r ˙
Therefore, to make M z = 0 , we need r ˙ = 0 . Assuming the initial state r 0 = 0 , this means r ˙ = r = 0 for the whole trajectory. Based on this, we can modify the inverse dynamics analysis to solve for ψ ˙ and ψ ¨ as:
r = ψ ˙ + sin θ cos ϕ p cos ϕ cos θ = 0 ψ ˙ = sin θ cos ϕ p
r ˙ = cos ϕ cos θ ψ ¨ + sin θ cos ϕ p ˙ + p cos θ cos ϕ θ ˙ + sin θ sin ϕ ϕ ˙ cos 2 ϕ r cos θ sin ϕ ϕ ˙ sin θ cos ϕ θ ˙ cos 2 ϕ = 0 ψ ¨ = sin θ cos ϕ p ˙ p cos θ cos ϕ θ ˙ + sin θ sin ϕ ϕ ˙ cos 2 ϕ
Having ψ ˙ and ψ ¨ at each moment in the trajectory, we can integrate ψ iteratively throughout the inverse dynamics analysis by ψ n + 1 = ψ n + ψ ˙ Δ t + 1 2 ψ ¨ Δ t 2 . Figure 19 and Figure 20 show a simple four-waypoint aggressive trajectory and the corresponding optimal yaw trajectory obtained by this method.
We converted this optimal yaw trajectory into a polynomial form and used the geometric controller to track this aggressive polynomial trajectory. Figure 21 shows that, by tracking this optimal yaw trajectory, the M z control input will be indeed close to zero.
Figure 22 shows the comparison of rotor thrust used to track the aggressive trajectory between constant zero yaw and the optimal yaw trajectory. The additional yaw control effort M z in the constant zero yaw case pulls the F 1 and F 3 rotor forces away from the F 2 and F 4 rotor forces and thereby increases the maximum rotor force needed. From the result, we can see that, by tracking the optimal yaw trajectory, the maximum rotor force is reduced significantly by 17.5% (from 4 N to 3.3 N).
Remark 1.
With the optimization and analysis framework developed, we review the choice of minimum snap trajectory. For the same four-waypoint scenario and segment times, we compare the optimal trajectory and the maximum rotor force required with the cost of different orders of trajectory derivatives to minimize. Figure 23 shows the geometry of the trajectories. As shown in Figure 24, the trajectory that minimizes the cost on the fourth derivative has the minimum force required. Therefore, we conclude that minimum snap is indeed the optimal choice for this approach of multi-copter aggressive trajectory optimization.

8. Conclusions

With the framework developed in this research, we can find the optimal polynomial trajectory and the corresponding segment time allocation for a given quadcopter model, a specified scenario, and a desired aggressiveness. We can also tell if the given scenario is beyond the capability of the given quadcopter by examining whether the solution exists. Furthermore, we can use the algorithm to evaluate the minimum required rotor thrust to achieve the scenario. Instead of tracking the maximum force, we can also use the algorithm to track the total time. That is, given a desired flight time, we can find the optimal trajectory that has the minimum required rotor thrust. Though the quadcopter model is used in this paper, the inverse dynamics analysis for control input (collective force and three-axis moment) applies to generic multi-copters. Therefore, the method developed in this research can be applied to multi-copters with minor modifications to the vehicle control allocation model. The geometric controller used in this research showed the capability of tracking the optimized aggressive trajectory. The tracking result verified the feasibility of the optimized trajectory and the credibility of the inverse dynamics analysis. In addition, the yaw trajectory optimization method is capable of improving the aggressive performance in the case of no requirement on heading angle.

Author Contributions

Both K.S. and T.-L.L. were involved in the conceptualization and development of methodology. T.-L.L. performed the simulations and wrote the article. K.S. reviewed and edited the manuscript for publication. All analyses were jointly performed by T.-L.L. and K.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Office of Naval Research under grant number N00014-18-1-2215 and The National Science Foundation (S&AS) under grant number 1724248. The APC was funded by Kamesh Subbarao.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors acknowledge the support from the Office of Naval Research under grant number N00014-18-1-2215 and the National Science Foundation (S&AS) under grant number 1724248.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Tal, E.; Karaman, S. Accurate Tracking of Aggressive Quadrotor Trajectories Using Incremental Nonlinear Dynamic Inversion and Differential Flatness. IEEE Trans. Control. Syst. Technol. 2020, 29, 1203–1218. [Google Scholar] [CrossRef]
  2. Gao, F.; Wu, W.; Pan, J.; Zhou, B.; Shen, S. Optimal Time Allocation for Quadrotor Trajectory Generation. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, Madrid, Spain, 1–5 October 2018. [Google Scholar]
  3. Yu, G.; Cabecinhas, D.; Cunha, R.; Silvestre, C. Quadrotor trajectory generation and tracking for aggressive maneuvers with attitude constraints. IFAC-PapersOnLine 2019, 52, 55–60. [Google Scholar] [CrossRef]
  4. Richter, C.; Bry, A.; Roy, N. Polynomial Trajectory Planning for Quadrotor Flight. In Proceedings of the International Conference on Robotics and Automation, Karlsruhe, Germany, 6–10 May 2013. [Google Scholar]
  5. Mellinger, D.; Kumar, V. Minimum snap trajectory generation and control for quadrotors. In Proceedings of the International Conference on Robotics and Automation, Shanghai, China, 9–13 May 2011. [Google Scholar]
  6. Mellinger, D.W. Trajectory Generation and Control for Quadrotors. Publicly Access. Penn Diss. 2012, 547, 5–9. [Google Scholar]
  7. Richter, C.; Bry, A.; Roy, N. Polynomial Trajectory Planning for Aggressive Quadrotor Flight in Dense Indoor Environments. In Robotics Research; Inaba, M., Corke, P., Eds.; Springer International Publishing: Cham, Switzerland, 2016; pp. 649–666. [Google Scholar]
  8. Bertsekas, D.P. Nonlinear Programming; Athena Scientific: Belmont, MA, USA, 1999. [Google Scholar]
  9. Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, MA, USA, 2004. [Google Scholar]
  10. Quan, Q. Introduction to Multicopter Design and Control; Springer: Singapore, 2017. [Google Scholar]
  11. Lee, T.; Leok, M.; McClamroch, N. Geometric tracking control of a quadrotor uav on SE(3). In Proceedings of the 49th IEEE Conference on Decision and Control, Atlanta, GA, USA, 15–17 December 2010. [Google Scholar]
Figure 1. Reference frames and quadcopter forces/moments.
Figure 1. Reference frames and quadcopter forces/moments.
Aerospace 09 00281 g001
Figure 2. Methodology of the aggressive trajectory optimization.
Figure 2. Methodology of the aggressive trajectory optimization.
Aerospace 09 00281 g002
Figure 3. Relevant frames for quadcopter and window constrains.
Figure 3. Relevant frames for quadcopter and window constrains.
Aerospace 09 00281 g003
Figure 4. Comparison of trajectories optimized with different c T .
Figure 4. Comparison of trajectories optimized with different c T .
Aerospace 09 00281 g004
Figure 5. Comparison of maximum rotor force and total time for trajectories optimized with different c T .
Figure 5. Comparison of maximum rotor force and total time for trajectories optimized with different c T .
Aerospace 09 00281 g005
Figure 6. Plot of the test scenario and the optimal trajectory.
Figure 6. Plot of the test scenario and the optimal trajectory.
Aerospace 09 00281 g006
Figure 7. Vehicle Euler angles along the trajectory flight.
Figure 7. Vehicle Euler angles along the trajectory flight.
Aerospace 09 00281 g007
Figure 8. Vehicle body angular rate along the trajectory flight.
Figure 8. Vehicle body angular rate along the trajectory flight.
Aerospace 09 00281 g008
Figure 9. Vehicle body angular acceleration along the trajectory flight.
Figure 9. Vehicle body angular acceleration along the trajectory flight.
Aerospace 09 00281 g009
Figure 10. Control input along the trajectory flight.
Figure 10. Control input along the trajectory flight.
Aerospace 09 00281 g010
Figure 11. Rotor thrust along the trajectory flight.
Figure 11. Rotor thrust along the trajectory flight.
Aerospace 09 00281 g011
Figure 12. Plot of the quadcopter along the trajectory flight.
Figure 12. Plot of the quadcopter along the trajectory flight.
Aerospace 09 00281 g012
Figure 13. Trajectory tracking result.
Figure 13. Trajectory tracking result.
Aerospace 09 00281 g013
Figure 14. Vehicle Euler angle comparison.
Figure 14. Vehicle Euler angle comparison.
Aerospace 09 00281 g014
Figure 15. Vehicle linear acceleration tracking.
Figure 15. Vehicle linear acceleration tracking.
Aerospace 09 00281 g015
Figure 16. Vehicle angular acceleration comparison.
Figure 16. Vehicle angular acceleration comparison.
Aerospace 09 00281 g016
Figure 17. Control input comparison.
Figure 17. Control input comparison.
Aerospace 09 00281 g017
Figure 18. Rotor force comparison.
Figure 18. Rotor force comparison.
Aerospace 09 00281 g018
Figure 19. Simple 4-waypoint aggressive trajectory.
Figure 19. Simple 4-waypoint aggressive trajectory.
Aerospace 09 00281 g019
Figure 20. Optimal yaw trajectory for the trajectory.
Figure 20. Optimal yaw trajectory for the trajectory.
Aerospace 09 00281 g020
Figure 21. Simulation result of control input with geometric controller.
Figure 21. Simulation result of control input with geometric controller.
Aerospace 09 00281 g021
Figure 22. Rotor force comparison.
Figure 22. Rotor force comparison.
Aerospace 09 00281 g022
Figure 23. Comparison of trajectories optimized with cost on different order of derivative.
Figure 23. Comparison of trajectories optimized with cost on different order of derivative.
Aerospace 09 00281 g023
Figure 24. Maximum rotor force required for trajectories optimized with cost on different order of derivative.
Figure 24. Maximum rotor force required for trajectories optimized with cost on different order of derivative.
Aerospace 09 00281 g024
Table 1. Quadcopter parameters.
Table 1. Quadcopter parameters.
m1.023 kgg 9.81 m / s 2
L0.2223 m I x x 0.0095 kg · m 2
k F 1.4865 · 10 7 N / RPM 2 I y y 0.0095 kg · m 2
k M 2.9250 · 10 9 N · m / RPM 2 I z z 0.0186 kg · m 2
Table 2. Waypoint settings in the scenario.
Table 2. Waypoint settings in the scenario.
WPT Setting x ( m ) y ( m ) z ( m ) ψ
Waypoint 10200
Waypoint 2120-
Waypoint 3100.5-
Waypoint 4000.5 180
Table 3. Window (WDW) settings in the scenario.
Table 3. Window (WDW) settings in the scenario.
WDW Setting ϕ θ ψ Location
Window 1 0 15 0 Waypoint 2
Window 2 30 0 20 Waypoint 3
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, T.-L.; Subbarao, K. Optimal Aggressive Constrained Trajectory Synthesis and Control for Multi-Copters. Aerospace 2022, 9, 281. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9060281

AMA Style

Liu T-L, Subbarao K. Optimal Aggressive Constrained Trajectory Synthesis and Control for Multi-Copters. Aerospace. 2022; 9(6):281. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9060281

Chicago/Turabian Style

Liu, Tsung-Liang, and Kamesh Subbarao. 2022. "Optimal Aggressive Constrained Trajectory Synthesis and Control for Multi-Copters" Aerospace 9, no. 6: 281. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9060281

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