In the present work, the tail-sitter UAV (SWAN K1 PRO) is utilized to validate our proposed altitude-hold transition method, as shown in Figure 1. The SWAN K1 PRO is developed by HEQ UAV Technical Company in Shenzhen, China. The vehicle’s mass and wingspan are 1.3328 kg and 1.085 m, respectively. It lacks control surfaces, relying solely on the four propellers on the fuselage to generate all torques and forces. The inertial frame O x y z is defined as North-East-Down (NED). The body frame O b x b y b z b is defined as Front-Right-Down (FRD). O is the origin of the world coordinate frame, while O b is the center of gravity of the vehicle.

The states of the tail-sitter UAV are denoted by x full = p , v , R , ω . p ∈ R 3 and v ∈ R 3 are the vehicle position and velocity in the world-fixed reference frame, respectively; ω ∈ R 3 is the angular velocity in the body frame. R ∈ S O ( 3 ) denotes the rotation matrix from the body frame to the inertial frame; The inputs for the tail-sitter in terms of force and torques are u full = f , τ , where f and τ ∈ R 3 denote the thrust and control moment vector produced by four propellers, respectively.

Within these definitions, the translational and rotational dynamics of the quadrotor tail-sitter UAV can be described by the following equations: p ̇ = v (1) v ̇ = g + 1 m f R e 1 + R f a (2) R ̇ = R ω (3) J ω ̇ = τ + M a − ω × J ω (4)

Where g = ( 0,0,9.8 ) T is the gravitational acceleration; and m is the vehicle mass. e 1 = ( 1,0,0 ) T , e 2 = ( 0,1,0 ) T and e 3 = ( 0,0,1 ) T are the unit vectors. f a ∈ R 3 and M a ∈ R 3 are the aerodynamic force and moment in the body frame, respectively, and will be introduced in Section Modeling of Aerodynamics. J ∈ R 3 × 3 is the inertia tensor matrix of the vehicle about the body frame with the center of gravity as the origin. The notation ⋅ converts a 3D vector into an antisymmetric matrix. In this study, the direction of thrust is assumed to align with the x b . In a situation where propellers are installed at a fixed angle, it only requires a simple transformation using a constant matrix.

The present study uses an analytical aerodynamic model to determine the aerodynamic forces and moments of an aircraft. Generally, these forces and moments can be described as a set of functions depending on the AoA α and the sideslip angle β .

The aerodynamic force f a is modeled in the body frame as follows: f a = − cos ⁡ α 0 sin ⁡ α 0 1 0 − sin ⁡ α 0 − cos ⁡ α D Y L (5) where L , D , and Y are the lift, drag, and side force produced by the fuselage and wings, respectively. L , M , and N are the rolling, pitching, and yawing moments, respectively, along the body axes x b , y b , and z b . Hence, the aerodynamic moment vector M a can be modeled as: M a = L , M , N T (6)

According to Etkin and Reid [16], the aerodynamic forces and moments can be parameterized as: L = 1 2 ρ V 2 S C L ,   D = 1 2 ρ V 2 S C D ,   Y = 1 2 ρ V 2 S C Y L = 1 2 ρ V 2 S b C l ,   M = 1 2 ρ V 2 S c ̄ C m ,   N = 1 2 ρ V 2 S b C n (7) where V = v a is the magnitude of the air speed, ρ is the density of the air, S is the reference area of the wing, c ̄ is the mean aerodynamic chord, and b is the span of airplane. It should be noted that v a = v − w , where w is the wind velocity. C L , C D , and C Y are the lift coefficients, drag coefficients, and side force coefficients of the vehicle, respectively, while C l , C m , and C n are the rolling, pitching, and yawing moment coefficients, respectively. Each set of non-dimensional numbers is exclusively linked to the aircraft’s shape and attitude.

The total aerodynamic force f a in Equation 5 can be rewritten as: f a = 1 2 ρ V 2 S C α , β (8) where C α , β = C x α , β C y α , β C z α , β T (9) C x α , β = − C D α , β cos ⁡ α + C L α , β sin ⁡ α (10) C y α , β = C Y α , β (11) C z α , β = − C D α , β sin ⁡ α − C L α , β cos ⁡ α (12)

Based on the air velocity with respect to the body frame v a B = R T v a = ( v a x B , v a y B , v a z B ) T , the AoA α and the sideslip angle β can be calculated as: α = arctan v a z B v a x B (13) β = arcsin v a y B V (14)

Since the airframe is symmetric with respect to the x b O b z b plane, it can be concluded that: C L α , β = C L α , − β , ∀ α , β (15) C D α , β = C D α , − β , ∀ α , β (16) C Y α , β = − C Y α , − β , ∀ α , β (17)

Computational fluid dynamics (CFD) are used to obtain the aerodynamic parameters. Details of the CFD method and simulation results can be found in our previous work [17]. Figure 2 shows an example of the surface pressure on the tail-sitter.

The essence of the altitude-hold transition can be regarded as an optimal boundary value problem (OBVP) that considers only the initial and final states, which can be expressed as follows: min ∫ 0 T p s t T p s t d t (18) s . t . t ∈ 0 , T (19) p s − 1 0 = s 0 , p s − 1 T = s f (20) p ̇ t < V max , p ̈ t < a max (21)

Where T is used to denote the temporal duration of the trajectory, and p ( s ) ( t ) means that we minimize the s -th derivative of position p . V max and a max are the maximum speed constraint and maximum acceleration constraint, respectively. s 0 , s f ∈ R s × 3 are the constant boundary conditions. p [ s − 1 ] ∈ R s × 3 is defined by: p s − 1 = p , p ̇ , … , p s − 1 T (22)

In practice, minimizing snap ( s = 4 ) optimization roughly corresponds to reducing the required control moment [18]. This increases the likelihood that the control input limits are satisfied and that the trajectory is feasible. Thus, the minimum snap trajectory is utilized in this work. Traditionally, an unconstrained OBVP for trajectory generation can be solved using Pontryagin’s Maximum Principle (PMP). The MINCO (minimum control) trajectory class [19] utilizes optimality conditions to directly calculate an OBVP. If the transfer time between the two points is prescribed a priori, the minimum control effort trajectory can be derived in closed form from the structure of the problem, without solving an optimal control problem. Specifically, for a system controlled by the s -th derivative of position, the trajectory is a polynomial of degree 2 s − 1 with respect to time, where the expression for a 3-D trajectory is denoted by: p t = c T ρ t , t ∈ 0 , T (23) where c ∈ R 2 s × 3 are the coefficients of ( 2 s − 1 ) th-degree polynomials, ρ ( x ) = ( 1 , x , x 2 , … , x ( 2 s − 1 ) ) T . To acquire the matrix c , we establish a linear system: A c = b (24) where A ∈ R 2 s × 2 s and b ∈ R 2 s × 3 are: A = H 0 0 0 G M (25) b = s 0 T , s f T T (26) H 0 , G M ∈ R s × 2 s are defined as: H 0 = ρ 0 , ρ ̇ 0 , … , ρ s − 1 0 T (27) G M = ρ T , ρ ̇ T , … , ρ s − 1 T T (28) The matrix A is nonsingular and banded for any T > 0 . Using the MINCO trajectory class Equations 18–21, can be recast as an unconstrained optimization problem with respect to the decision variable T . Substituting the optimal T ∗ into Equation 24 and solving then yields a minimum-snap trajectory.

Leveraging the differential flatness property, several recent studies [20, 21] for multirotor UAVs have explicitly incorporated attitude, angular-velocity, and control-input constraints already at the trajectory generation phase. The underlying reason is that, for multirotors, both the system states and the control inputs can be expressed as closed-form functions of position and a finite number of its time derivatives, thereby enabling the trajectory optimization to be cast as a second-order cone program (SOCP) and solved efficiently. In contrast, in the differential flatness mapping of a tail-sitter UAV, the attitude must be obtained by numerically solving nonlinear algebraic equations, making it difficult to derive a closed-form expression in terms of the decision variable T . Consequently, attitude cannot be included directly in the optimization as a function of T , nor can angular-velocity and control-input constraints be explicitly handled as convex constraints at the planning stage, which hinders the formulation of the problem as an SOCP. Nevertheless, this does not imply a deficiency in the feasibility of the proposed trajectory-generation method: we implicitly bound the states by imposing a maximum-speed constraint and a maximum-acceleration constraint, and we implicitly constrain variations in the control inputs by minimizing the second time derivative of acceleration (snap). A substantial body of prior work [9, 22, 23] has shown that enforcing constraints on velocity and its finite-order derivatives (equivalently, on higher-order derivatives of position) at the planning stage suffices, to a large extent, to ensure dynamic feasibility. In addition, during the trajectory tracking phase, control-input constraints will be explicitly enforced (see Section Global Controller for the Tail-Sitter).

The system Equations 1–4 is a cascaded structure; that is, the input torque, τ , only affects the angular velocity. The attitude of the UAV is then calculated through the angular velocity, and the current speed and position can be determined through the attitude of the UAV. Due to this cascaded dynamics, the angular velocity dynamics Equation 4 can be tracked separately. In this work, the angular velocity of the tail-sitter is tracked by three decoupled proportional–integral–derivative (PID) controllers in the autopilot. The Coriolis term, ω × J ω , and the aerodynamic moment, M a , are both regarded as unknown disturbances and are compensated for in a feed-forward way. With a well-designed PID controller, the vehicle’s angular velocity can be achieved instantaneously. Thus, in terms of controller design, the tail-sitter’s states and inputs can be simplified as follows: x = p T v T R T T (29) u = f ω T T (30)

Based on Equations 29, 30, we utilize the same controller as in our previous work [17], which is derived from Model Predictive Control (MPC). MPC is a common and powerful trajectory tracking method in robotics. Generally, we can design MPC to track the trajectory of a model whose states are all in Euclidean space (such as the unicycle model). This is because all system states are flat vectors, and by calculating their error’s norm, we can quantify the tracking error. However, for the dynamic model shown in Equations 1–4, the tail-sitter’s states contain the rotation matrix R ∈ R 3 × 3 , which is not in Euclidean space but on a non-vector, curved manifold. This makes it difficult to quantify the rotation matrix error δ R . Inspired by [24], the rotation matrix R can be projected into Euclidean space via mapping: θ = L o g R ∈ R 3 (31) where L o g ( ⋅ ) is the logarithmic map proposed by [25], which is defined as: L o g S O 3 R = ϕ 2 ⁡ sin ⁡ ϕ R 32 − R 23 R 13 − R 31 R 21 − R 12 (32) where ϕ satisfies cos ( ϕ ) = 1 2 ( t r ( R ) − 1 ) .

Furthermore, δ R can be formulated as: δ R = R r ⊟ R = L o g R T R r (33) where R r is the rotation matrix on the reference trajectory, and R is the actual rotation matrix. The error of the rotation matrix is represented as a 3D vector, and calculating its norm can track the attitude of the vehicle.

Next, we construct the MPC objective function to ensure optimal performance of our control strategy. Benefiting from the differential flatness of quadrotor tail-sitter UAVs, the UAV’s reference 12-state can be determined from a time differentiable curve at any given moment. Therefore, according to Equations 29, 30, 33, the error of the states and inputs can be represented as: δ x = δ p T δ v T δ R T T ∈ R 9 δ u = δ f δ ω T T ∈ R 4 (34) δ p = p r − p ∈ R 3 (35) δ v = v r − v ∈ R 3 (36) δ R = R r ⊟ R = L o g R T R r ∈ R 3 (37) δ f = f r − f ∈ R (38) δ ω = ω r − ω ∈ R 3 (39)

Where the subscript r represents the value on the reference trajectory, and no subscript represents the actual value. Therefore, we can construct an MPC as an optimization problem: δ u k * = arg min δ u k ∑ k = 0 N − 1 x r k + 1 − x k + 1 Q k 2 + δ u k P k 2  s.t.  δ x k + 1 = I + Δ t F x k δ x k + Δ t F u k δ u k + Δ t F w k δ w k u min ≤ u r k − δ u k ≤ u max (40) where N is the receding horizon, Q k and P k are positive diagonal matrices, u min and u max denote the actual input constraints, F x , F u , and F w are the Jacobians of the error-state dynamic model δ x ̇ with respect to δ x , δ u and δ w , respectively. The δ w = w r − w , where w is the wind velocity in the environment, and w r is the wind velocity used in trajectory optimization. The objective function should be as small as possible so that the states of the trajectory are well tracked and the trajectory inputs are as smooth as possible. Finally, the optimal control of a trajectory is: u * = u r − δ u * (41)

To obtain F x , F u , and F w , we take the derivative of the state error with respect to time to obtain the error dynamic model: δ p ̇ = δ v (42) δ v ̇ = 1 m f r R r e 1 + R r f a r − f R e 1 + R f a (43) δ R ̇ = A − T δ R − R r T R ω + ω r (44)

Where f a r and f a are the aerodynamic forces on the reference trajectory and actual trajectory, respectively. A ( ⋅ ) : R 3 → R 3 × 3 represents a map [25]: A δ R = I + 1 − cos ‖ δ R ‖ ‖ δ R ‖ ⌊ δ R ⌋ ‖ δ R ‖ + 1 − sin ‖ δ R ‖ ‖ δ R ‖ ⌊ δ R ⌋ 2 ‖ δ R ‖ 2 (45)

See Supplementary Appendix SA for the proof of Equation 44.

Linearizing Equations 42–44: δ x ̇ = F x δ x + F u δ u + F w δ w (46) where F x ≈ 0 I 0 0 1 m R r ∂ f a r ∂ v a B r R r T ∂ δ v ̇ ∂ δ R 0 0 − ω − 1 2 δ ω + 1 2 K (47) F u ≈ 0 0 R r e 1 0 0 I + 1 2 δ R (48) F w ≈ 0 − 1 m R r ∂ f a r ∂ v a B r R r T 0 (49) ∂ δ v ̇ ∂ δ R ≈ R r − a T e 1 − 1 m f a + 1 m ∂ f a r ∂ v a B r R r T v a (50)

Where the matrix K and the proof of Section Global Controller for the Tail-Sitter are given in Supplementary Appendix SB.

In conclusion Equation 40, can be seen as an optimization problem with an inequality constraint that can be solved by the Powell-Hestenes-Rockafellar augmented Lagrangian method (PHR-ALM) [26]. The controller structure is shown in Figure 3.