# VTOL UAV Transition Maneuver Using Incremental Nonlinear Dynamic Inversion

Authors: Zhenchang Liu; Jie Guo; Mengting Li; Shengjing Tang; Xiao Wang
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

The paper seeks to study the control system design of a novel unmanned aerial vehicle (UAV). The UAV is capable of vertical takeoff and landing (VTOL), transition flight and cruising via the technique of direct force control. The incremental nonlinear dynamic inversion (INDI) approach is adopted for the 6-DOF nonlinear and nonaffine control of the UAV. Based on the INDI control law, a method of two-layer cascaded optimal control allocation is proposed to handle the redundant and coupled control variables. For the weight selection in optimal control allocation, a dynamic weight strategy is proposed. This strategy can adjust the weight of the objective function according to the flight states and mission requirements, thus determining the optimizing direction and ensuring the rationality of the allocation results. Simulation results indicate that the UAV can track the target trajectory accurately and exhibit continuous maneuverability in transition flight.


---

## Body

## 1. Introduction

UAVs have increasing applications including surveillance, communications, search and rescue operations, and other military tasks. Among different flight conditions of UAVs, the 0–1000 m height area in cities is one of the most significant applications, where complex terrain and significant gusts arising from atmospheric turbulence exist.This research studies a novel fixed-wing VTOL UAV with thrust vector engines, which can be applied in urban areas. The UAV can get rid of the restrictions imposed by takeoff and landing conditions and be accurately recovered by hover function; the UAV owns a larger combat range and higher flight speed by forward flight capability. A new concept of low-speed cruising is studied in this paper. For the most current VTOL aircrafts, the transition from hover to forward flight is short and stable. However, for low-speed cruising, the transition is prolonged as a normal flight state. The vehicle can maintain transition flight for long time cruising by adopting direct force control and has favorable maneuverability. This flight mode is appropriate for vehicles flying in low-altitude complex conditions. In the meantime, under the low speed of UAV and the inefficient aerodynamic surface control during the transition flight, the incorporated control strategy of vectoring nozzles and aerodynamic surface should be adopted to control the attitude. In view of this, the thrust vector direct force control is of critical significance for transition maneuver. However, large nonlinearities, redundancies, and coupling effects arise when this technique is adopted.In recent years, the research on VTOL UAVs is increasingly prosperous with the advancements in automatic control and the increasing popularity of UAV platforms [1]. The problem of transition maneuvers has been studied for different types of UAVs, including a fixed-wing aircraft equipped with a thrust vector engine and lift fan [1–3], tilt-rotor aircraft [4–7], tail-sitter aircraft [8, 9], ducted-fan VTOL aircraft [10, 11], and tilt-wing aircraft [12, 13]. A pilot auxiliary control system was designed by Francesco and Mattei [1] for a fixed-wing tiltrotor UAV, and the control logics for different flight states were specified and synthesized through adopting INDI method. The daisy-chaining logic was employed to handle the control redundancy. The autonomous transition control of two V/STOL aircrafts was studied by Xili et al. [2]. The nonlinear trajectory control strategies in longitudinal direction have been designed specifically for aircrafts in different types. The maneuver of a tilt quadcopter was researched by Ryll et al. [7]. The control design was basically dependent on the exact linearization of the motion equations, and the actuation redundancy was calculated through employing pseudoinverse matrices. The previous works generally decouple the longitudinal and lateral control, design control logics for hover, transition flight, and forward flight, respectively, and regard some redundant control variables as constants to solve the transition maneuver control problem. The UAV’s maneuver potential with vectored thrust cannot be fully utilized.The INDI method is adopted in this paper to control UAV’s position/attitude during transition maneuver. The INDI method, which originates from the nonlinear dynamic inversion (NDI), solves the incremental form of equations of motion and generates a control law substantially reducing the dependence on aerodynamic model and other vehicle models. INDI was firstly adopted to control UAV’s attitude control by Sieberling et al. [14]. The INDI attitude control law of a quadcopter was proposed, and the momentum of the propellers was incorporated in the controller by Smeur et al. [15]. They generalized this method based on the previous work to the outer control loop (position control loop) of a quadcopter under severe gust loads [16, 17]. Lu et al. [18] applied the INDI method to the fixed-wing aircraft trajectory controller. They compared the performance of the INDI controller with that of NDI approach and proved that the INDI method can reduce the model uncertainties effectively.An INDI control system is designed in this paper to address the 6-DOF nonlinear control of the UAV in transition flight, and the main contributions are listed below:(1)
Given the problems of strong nonlinearities and multiaxis coupling characteristics, a unified 6-DOF nonlinear control strategy is proposed to control position/attitude and there is no need to switch the control logic according to different flight states. The INDI method is introduced to address the model uncertainty and control coupling problem. Different from the work conducted by Lu et al. [18], the sideslip angle β is not assumed to be zero, and the vectored thrust in 3 directions of the body axis is considered
(2)
A two-layer cascaded optimal control allocation method is proposed to address the control redundancy based on the INDI control law. The first-layer optimal allocation is conducted to allocate the increment of flight attitude and vectored thrust in the translational dynamics control loop. The solution of the engine thrust, vectoring nozzle deflections, and aerodynamic surface deflections are calculated in the second-layer control allocation
(3)
A dynamic weight selection strategy is designed for the objective function of the two-layer cascaded optimal control allocation. In dynamic weight selection strategy, a weight generator and a weight regulator are designed, which calculate weight through an analytic hierarchy process (AHP) and adjust weight according to flight states and mission requirement, respectively. The dynamic weight selection strategy can allocate control variables according to mission requirement and ensure the optimal results to be feasibleThis paper is constructed as follows: the configuration and aerodynamics characteristics of the VTOL UAV researched in this paper are described in Section2. A 6-DOF mathematical model and an INDI control system of the vehicle are given in Section 3. The two-layer cascaded optimal control allocation method is presented in Section 4. Section 5 provides the simulations results. Eventually, conclusions are drawn in Section 6.


## 2. UAV Configuration and Aerodynamics

The UAV is designed in a tandem-wing plus lift-body configuration. This aerodynamic configuration can provide more lift under the limitation imposed by wing span, and it makes UAV applicable to fly in low-altitude complex flight condition. The power system consists of a lift fan in the front part of fuselage and two thrust vector engines in each side of the rear part of fuselage, as presented in Figures1(a) and 1(b). At the bottom of the lift fan, a control rudder is equipped with the shaft. The control rudder can be deflected left and right with 12 degrees and provides lateral vectored thrust Ty and yaw moments nT, as illustrated in Figures 1(c) and 1(d). The engine’s vectoring nozzle can be swung up with 15 degrees and down with 90 degrees. It provides 3-axis vectored thrust for a direct force control, as shown in Figures 1(e) and 1(f). Through the cooperation of the lift fan and thrust vector engine, the UAV can realize VTOL, transition flight, and cruise.Figure 1
UAV aerodynamic and power configurations.
(a)
Vectored thrust of power system
(b)
UAV power system
(c)
Vectored thrust of lift fan
(d)
Lift fan
(e)
Vectored thrust of engine
(f)
Thrust vector engineIn the case of no wind tunnel experiment, the UAV’s aerodynamic coefficients were obtained by the method of aerodynamic estimation and CFD calculation together. In aerodynamic coefficients, the static pitch moment derivativeCmα ranges from −0.04 to 0.056. When the angle of attack is from 0 to 18 degrees, Cmα>0, which indicates that the UAV is statically instable in pitch channel. The reason of instability design is that for VTOL UAV, the alignment of the mass center should be considered with power system distribution and the location of aerodynamic center, which is in the front of the mass center, facilitates the vehicle to pitch up in transition maneuver. The static yaw moment derivative Cnβ ranges from −0.017 to −0.03. The static stability in yaw channel is caused by the rear wing dihedral effect. The vehicle’s main dynamic moment derivatives are Clp=−0.0084, Cmq=−0.01739, and Cnr=−0.0011. They are all less than zero, and it indicates that the damping moments can reduce the angular velocities and keep the vehicle dynamics stable. The mapping relationship between the actuators, aerodynamic forces, and moments are usually nonlinear. In this paper, to meet the requirement of the control allocation, polynomial approximation is used to fit the aerodynamic data curve between the actuator deflections and aerodynamic coefficients. This will alleviate the need of large lookup tables and speed up the control allocation computational process. The basic data of UAV is presented in Table 1.Table 1
Prototype basic data.
Item
Data
Running takeoff speed
15 m/s
Cruise speed
20 m/s
Vehicle length
1 m
Rear wing span
1.57 m
Vehicle weight
5 kg
Maximum engine thrust (single)
2.6 kg
Maximum lift fan thrust
3.8 kg
Flight range
6 km


## 3. INDI Control System Design

This section presents the design process of the control system. The trajectory tracking controller, by adopting time-scaled method, is split into four control loops, i.e., translational kinematics (position) control loop, translational dynamics (flight path) control loop, rotational kinematics (attitude) control loop, and rotational dynamics (angular rate) control loop. Given the differences of the inner and outer loops in the time constants, the control laws for inner and outer loops can be designed independently [19, 20]. Considering that the model uncertainties caused by aerodynamic force and moments merely exist in translational dynamics and rotational dynamics control loops, these two control loops are designed by the INDI method while the other two loops are designed by the NDI method. The equations of the 6-DOF flight dynamic model and thrust vector engine model adopted to design the control system are described in this section.


### 3.1. Translational Kinematic Control Loop

The reference flight path vector is calculated in translational kinematic control loop based on the target trajectory. The position vector and flight path vector of the aircraft are defined as(1)x0=xyzT,x1=VχγT,where x, y, and z represent the positions of the aircraft in the North, East, and down directions, respectively. V denotes the total velocity of the aircraft, χ indicates the kinematic azimuth angle, and γ represents the flight path angle. The equation of translational kinematics is expressed as
(2)ẋ0=g0x1=VcosχcosγVsinχcosγ−Vsinγ.The desired derivatives of position vector can be designed by a classical linear controller according to the reference trajectory and the feedback of the vehicle’s position.(3)ẋ0des=ẋ0ref+K0x0−x0ref.In (3), K0=KxKyKzT are the gains of the position linear controller. For a better explanation, the definitions are given below: the variables with superscript “des” denote the desired commands generated by the linear controller, and the variables with superscript “ref” denote the reference commands given for the controller to follow. The translational kinematic control loop contains no model uncertainties, and the NDI method is adopted to calculate the control input. The reference flight path command for the next control loop is calculated in
(4)x1ref=−Żdessinarctan−Żdes/ẎdessinarctanẎdes/ẊdesarctanẎdesẊdesarctan−ŻdesẎdessinarctanẎdesẊdes.


### 3.2. Translational Dynamic Control Loop

#### 3.2.1. INDI Control Law Design

In translational dynamic control loop, the vehicle’s flight path statusx1 is controlled by aerodynamic force and vectored thrust together to track the reference command x1ref obtained by the last control loop. Define the control variables as
(5)x2=αβμTxTyTzT,x2a=αβμT,x2t=TxTyTzT.In this loop, there are totally 6 control inputs. Among them,x2a represent the aerodynamic control variables, where μ is the bank angle, and α and β are the angle of attack and angle of sideslip, respectively. x2t represent the vectored thrusts, and Tx, Ty, and Tz are the components of the thrust in the body axis. The equation of translational dynamics is expressed as
(6)ẋ1=f1x1+g1x1,x2,(7)f1x1=1m0001mVcosγ000−1mVLkg00mg,(8)g1x1,x2=1m0001mVcosγ000−1mVLkbα,β,μx2t+LkaμFA.In (7) and (8), Lkg is the direction cosine matrix from earth axis to flight path axis, Lkbα,β,μ is the direction cosine matrix from body axis to flight path axis, and Lkaμ is the direction cosine matrix from wind axis to flight path axis. FA=−DαCβ−LαT is the aerodynamic force in which D represents drag, C represents aerodynamic side force, and L represents lift. FA and Lkbα,β,μ and Lkaμ are assumed as the functions of x2a. The multiplication of direction cosine matrixes causes severe coupling in control variables. Furthermore, in translational dynamic control loop, the error caused by the estimation of FA will bring uncertainties in the control system. Therefore, the INDI method is adopted to address the control problem. To rewrite (6) in incremental form, it is defined that x1n and x2n represent the values of x1 and x2 in the next time step, and their relationships are described below:
(9)x1n=x1+Δx1,x2n=x2+Δx2.Applying Taylor expansion toẋ1n at x1,x2 and with higher-order terms neglected, the result can be expressed as
(10)ẋ1n=ẋ1+∂f1x1∂x1x1n−x1+∂g1x1,x2∂x1x1n−x1+∂g1x1,x2∂x2x2n−x2.In (10), the second and third terms partial to x1 are assumed much smaller than the forth term, partial to x2. This commonly arises from the principle of time scale separation [20]. For simplification, it is approximated that
(11)∂g1x1,x2∂x2x2n−x2=g1x1,x2+Δx2−g1x1,x2Δx2x2n−x2=g1x1,x2+Δx2−g1x1,x2.Replaceẋ1n with ideal derivatives of flight path vector ẋ1des, which is calculated by the linear controller. The incremental control equation can be denoted as
(12)ẋ1des−ẋ1=g1x1,x2+Δx2−g1x1,x2.To simplify the computation process, it is assumed that(13)Δx2iΔx2j=0,where x2i and x2j represent the elements in vector x2. Accordingly, the equation can be rewritten into the affine-in-control form
(14)ẋ1des−ẋ1=g1a00g1tΔx2aΔx2t,where g1a and g1t are the 3×3 matrixes, representing the aerodynamic control matrix and thrust vector control matrix. The detail information of g1a and g1t is shown in Appendix. Compared with a traditional aircraft in incremental approach, the thrust vector control matrix g1t is additional. The elements in g1t are polynomials of control variables in x2a. Thus, the control system is transferred into a linear and time variant system, and the increment of control variables is decoupled. Then the ideal incremental virtual command ∆x2 can be calculated by
(15)Δx2aΔx2t=g1a00g1t‡ẋ1des−ẋ1.In (15), the superscript “‡” represents the dynamic weight pseudoinverse method, which will be introduced in the next section. The reference command for the next control loop is obtained by
(16)x2ref=x2ax2t+Δx2aΔx2t.


#### 3.2.2. Flight Path Vector Derivative Acquisition

In (12), (14), and (15), ẋ1 can be derived by V, χ, and γ, which are measured by onboard sensors. With the ẋ1des−ẋ1 introduced, the term f1x1 is cancelled, which is the reason why INDI is referred to a sensor-based approach. This approach transfers the dependence on the model accuracy into the dependence on the sensor accuracy. In most cases, the signals were measured by sensor contain noise, and the differentiation of noisy signal amplifies the noise. To make results accurate, a filter can be adopted to abate the noise in sensor data. In this paper, a second-order filter is employed. As stated in literatures [16, 18, 21], the washout filter can be expressed in Laplace domain as
(17)Hs=ωn2s2+2ζωns+ωn2.However, the filter leads to a delay which should be compensated. In the Taylor expansion shown in (10), ẋ1, x1, and x2 should be from the same moment. In this regard, a second-order filter is also applied for x2 to counteract the impact caused by time delay in x1 and ẋ1, and (12) can be rewritten in a time-synchronized form
(18)ẋ1des−ẋ1f=g1x1f,x2f+Δx2−g1x1f,x2f,where subscript “f” represents the filtered variable. Alternatively, other methods can also deal with the measurement delay problem, such as predictive filtering proposed by Sieberling et al. [14]. However, the prediction requires additional modeling and cannot predict disturbances. The final reference command for x2ref can be given by
(19)x2ref=x2f+Δx2.


#### 3.2.3. Uncertainty Analysis

The model uncertainty stemmed from aerodynamic force in translational dynamic control loop is analyzed in this subsection. The change of aerodynamic coefficients is assumed to be primarily caused by the angle of attack (α) and sideslip angle (β) for simplification. The aerodynamic coefficient can be expressed as
(20)CD=C¯D0+CDαα,CC=C¯C0+CCββ,CL=C¯L0+CLαα.In (20), CD, CC, and CL represent the total drag coefficient, side force coefficient, and lift coefficient, respectively; C¯D0 represents the sum of drag coefficients without the part contributed by α; C¯C0 represents the sum of side force coefficients without the part contributed by β; C¯L0 represents the sum of lift coefficients without the part contributed by α; and CDα, CCβ, and CLα represent derivative coefficients about α and β. These coefficients vary with the flight condition and can be assumed as constants in the calculation of each control period. The methods to obtain these parameters involve mathematic estimation and CFD analysis, as described in Section 2. The errors between the estimate coefficients and accurate coefficients make the model uncertain and assume that the uncertainty is mainly caused by C¯D0, C¯C0, and C¯L0. As the work by Lu et al. [18] indicates, the kinematic roll angle (μ) can be calculated directly under the assumption of β=0, which largely simplifies the control equation. Additionally, coefficients C¯D0, C¯C0, and C¯L0 can be eliminated by adopting INDI, which evidently reduces model uncertainty. However, the UAV in this paper with direct force control owns great lateral maneuverability. Therefore, β cannot be approximated as zero, and μ is calculated through solving the control equation. As a result, coefficients C¯D0, C¯C0, and C¯L0 arise in the third column of g1a and cannot be eliminated. Nevertheless, no uncertainty exists in the thrust vector control matrix g1t.The application of INDI in the translational dynamics control loop has two primary advantages. First, INDI control law restrains the model uncertainty caused by aerodynamic force in control matrixg1a and isolates its impact on direct force incremental control. Second, INDI can decouple the control variables in g1x1,x2 and transfer the control equation into a linear form. Accordingly, the aerodynamic force and vectored thrust control allocation are simplified, the complex numerical solution of nonlinear coupling control is avoided, and computational load of onboard control system is reduced.


### 3.3. Rotational Kinematic Control Loop

The objective in this control loop is to trackx2aref given in the translational dynamics control loop. The angular rate vector of the vehicle is defined as
(21)x3=pqrT.In (21), p, q, and r represent the roll, pitch, and yaw rates in body axis, respectively. No model uncertainty exists in this control loop, and the control law based on the standard NDI approach is established below:
(22)ẋ2a=f2x2a+g2x2ax3,f2x2a=1cosβ−cosγ−sinμcosγ−sinμcosβcosμcosγcosβcosμsinβsinγcosβ+sinμcosγsinβγ̇χ̇,g2x2a=1cosβ−cosαsinβcosβ−sinαsinβsinαcosβ0−cosαcosβcosα0sinα.It is noteworthy thatγ̇, χ̇, and μ are not directly measured onboard in this control system. They are calculated by (23) derived from (6).
(23)γ̇=cosαsinβsinμ+sinαcosμAx−cosβsinμAy+sinαsinβsinμ−cosαcosμAzV,χ̇=−cosαsinβcosμ+sinαsinμAx+cosβcosμAy+−sinαsinβcosμ−cosαsinμAzVcosγ,where Ax, Ay, and Az are directly measured by the accelerometers of the aircraft. The reason why γ̇ and χ̇ are not acquired by the same method adopted in translational dynamic control loop is that the filter will delay the measurement and cause errors from accurate results. The reference angular rates can be obtained by
(24)x3ref=g2x2a−1ẋ2ades−f2x2a,where ẋ2ades is the desired command and designed by linear controller, which is similar to the one adopted to design ẋ1des.


### 3.4. Rotational Dynamic Control Loop

x4 is defined as the control moments acting on the vehicle, and it can be denoted as
(25)x4=x4s+x4t=lcmcncT,where x4s denote the moments generated by aerodynamic control surface; x4t indicate the moments generated by vectored thrust; and lc, mc, and nc represent the rolling, pitching, and yawing control moments, respectively. The dynamics of the angular rates of the vehicle can be expressed into the following affine-in-control form
(26)ẋ3=f3x3+J−1x4,f3x3=J−1Ma−x3×Jx3,where J denotes the inertia matrix and Ma refers to the aerodynamic moments generated by derivatives unrelated to control surface deflections.
(27)Ma=q¯SC¯l0C¯m0C¯n0.Like aerodynamic force, the aerodynamic coefficients can also be denoted as(28)Cl=C¯l0+Clδaδa+Clδrδr,Cm=C¯m0+Cmδeδe,Cn=C¯n0+Cnδaδa+Cnδrδr.In (27) and (28), Cl, Cm, and Cn represent the total 3-axis aerodynamic moment coefficients, respectively; C¯l0 represents the sum of aerodynamic rolling moment coefficients without the part contributed by δa and δr; C¯m0 represents the sum of aerodynamic pitching moment coefficients without the part contributed by δe; C¯n0 represents the sum of aerodynamic yawing moment without the part contributed by δa and δr; and Clδa, Clδr, Cmδe, Cnδa, and Cnδr represent derivative aerodynamic moment coefficients about δa, δe, and δr.The aerodynamic control moments are denoted as(29)x4s=q¯Sb1CASδAS,CAS=Clδa0Clδr0Cmδe0Cnδa0Cnδr,b1=diagb,c¯,b,where q¯ denotes the dynamic pressure, S denotes the reference wing area, b denotes the wing span, and c¯ denotes the mean aerodynamic chord. δAS=δaδeδrT are control surface deflections. The control moment coefficients CAS are assumed to be accurate in aerodynamic estimation, and the uncertainties are considered to be caused principally by C¯l0, C¯m0, and C¯n0. Rewriting (26) into the incremental form, as presented in (30), the term f3x3 can be cleared up.
(30)Δx4=Jẋ3des−ẋ3f.In (30), ẋ3f represents the derivative of x3, which is the filtered signal collected by the sensor. In the absence of f3x3, the uncertainty caused by Ma is eliminated and the nonlinear cross couplings of the angular rate term is also cancelled. The reference control moments can be calculated as
(31)x4ref=x4f+Δx4.After the calculation of four control loops,∆x2t and ∆x4 will be output as virtual commands to the second-layer control allocation. The control allocation of the engine thrust, vectoring nozzles, and aerodynamic control surface deflections will be introduced in the next section. The block diagram of the control system is illustrated in Figure 2.Figure 2
Conceptual block diagram of the control system.


## 3.1. Translational Kinematic Control Loop

The reference flight path vector is calculated in translational kinematic control loop based on the target trajectory. The position vector and flight path vector of the aircraft are defined as(1)x0=xyzT,x1=VχγT,where x, y, and z represent the positions of the aircraft in the North, East, and down directions, respectively. V denotes the total velocity of the aircraft, χ indicates the kinematic azimuth angle, and γ represents the flight path angle. The equation of translational kinematics is expressed as
(2)ẋ0=g0x1=VcosχcosγVsinχcosγ−Vsinγ.The desired derivatives of position vector can be designed by a classical linear controller according to the reference trajectory and the feedback of the vehicle’s position.(3)ẋ0des=ẋ0ref+K0x0−x0ref.In (3), K0=KxKyKzT are the gains of the position linear controller. For a better explanation, the definitions are given below: the variables with superscript “des” denote the desired commands generated by the linear controller, and the variables with superscript “ref” denote the reference commands given for the controller to follow. The translational kinematic control loop contains no model uncertainties, and the NDI method is adopted to calculate the control input. The reference flight path command for the next control loop is calculated in
(4)x1ref=−Żdessinarctan−Żdes/ẎdessinarctanẎdes/ẊdesarctanẎdesẊdesarctan−ŻdesẎdessinarctanẎdesẊdes.


## 3.2. Translational Dynamic Control Loop

### 3.2.1. INDI Control Law Design

In translational dynamic control loop, the vehicle’s flight path statusx1 is controlled by aerodynamic force and vectored thrust together to track the reference command x1ref obtained by the last control loop. Define the control variables as
(5)x2=αβμTxTyTzT,x2a=αβμT,x2t=TxTyTzT.In this loop, there are totally 6 control inputs. Among them,x2a represent the aerodynamic control variables, where μ is the bank angle, and α and β are the angle of attack and angle of sideslip, respectively. x2t represent the vectored thrusts, and Tx, Ty, and Tz are the components of the thrust in the body axis. The equation of translational dynamics is expressed as
(6)ẋ1=f1x1+g1x1,x2,(7)f1x1=1m0001mVcosγ000−1mVLkg00mg,(8)g1x1,x2=1m0001mVcosγ000−1mVLkbα,β,μx2t+LkaμFA.In (7) and (8), Lkg is the direction cosine matrix from earth axis to flight path axis, Lkbα,β,μ is the direction cosine matrix from body axis to flight path axis, and Lkaμ is the direction cosine matrix from wind axis to flight path axis. FA=−DαCβ−LαT is the aerodynamic force in which D represents drag, C represents aerodynamic side force, and L represents lift. FA and Lkbα,β,μ and Lkaμ are assumed as the functions of x2a. The multiplication of direction cosine matrixes causes severe coupling in control variables. Furthermore, in translational dynamic control loop, the error caused by the estimation of FA will bring uncertainties in the control system. Therefore, the INDI method is adopted to address the control problem. To rewrite (6) in incremental form, it is defined that x1n and x2n represent the values of x1 and x2 in the next time step, and their relationships are described below:
(9)x1n=x1+Δx1,x2n=x2+Δx2.Applying Taylor expansion toẋ1n at x1,x2 and with higher-order terms neglected, the result can be expressed as
(10)ẋ1n=ẋ1+∂f1x1∂x1x1n−x1+∂g1x1,x2∂x1x1n−x1+∂g1x1,x2∂x2x2n−x2.In (10), the second and third terms partial to x1 are assumed much smaller than the forth term, partial to x2. This commonly arises from the principle of time scale separation [20]. For simplification, it is approximated that
(11)∂g1x1,x2∂x2x2n−x2=g1x1,x2+Δx2−g1x1,x2Δx2x2n−x2=g1x1,x2+Δx2−g1x1,x2.Replaceẋ1n with ideal derivatives of flight path vector ẋ1des, which is calculated by the linear controller. The incremental control equation can be denoted as
(12)ẋ1des−ẋ1=g1x1,x2+Δx2−g1x1,x2.To simplify the computation process, it is assumed that(13)Δx2iΔx2j=0,where x2i and x2j represent the elements in vector x2. Accordingly, the equation can be rewritten into the affine-in-control form
(14)ẋ1des−ẋ1=g1a00g1tΔx2aΔx2t,where g1a and g1t are the 3×3 matrixes, representing the aerodynamic control matrix and thrust vector control matrix. The detail information of g1a and g1t is shown in Appendix. Compared with a traditional aircraft in incremental approach, the thrust vector control matrix g1t is additional. The elements in g1t are polynomials of control variables in x2a. Thus, the control system is transferred into a linear and time variant system, and the increment of control variables is decoupled. Then the ideal incremental virtual command ∆x2 can be calculated by
(15)Δx2aΔx2t=g1a00g1t‡ẋ1des−ẋ1.In (15), the superscript “‡” represents the dynamic weight pseudoinverse method, which will be introduced in the next section. The reference command for the next control loop is obtained by
(16)x2ref=x2ax2t+Δx2aΔx2t.


### 3.2.2. Flight Path Vector Derivative Acquisition

In (12), (14), and (15), ẋ1 can be derived by V, χ, and γ, which are measured by onboard sensors. With the ẋ1des−ẋ1 introduced, the term f1x1 is cancelled, which is the reason why INDI is referred to a sensor-based approach. This approach transfers the dependence on the model accuracy into the dependence on the sensor accuracy. In most cases, the signals were measured by sensor contain noise, and the differentiation of noisy signal amplifies the noise. To make results accurate, a filter can be adopted to abate the noise in sensor data. In this paper, a second-order filter is employed. As stated in literatures [16, 18, 21], the washout filter can be expressed in Laplace domain as
(17)Hs=ωn2s2+2ζωns+ωn2.However, the filter leads to a delay which should be compensated. In the Taylor expansion shown in (10), ẋ1, x1, and x2 should be from the same moment. In this regard, a second-order filter is also applied for x2 to counteract the impact caused by time delay in x1 and ẋ1, and (12) can be rewritten in a time-synchronized form
(18)ẋ1des−ẋ1f=g1x1f,x2f+Δx2−g1x1f,x2f,where subscript “f” represents the filtered variable. Alternatively, other methods can also deal with the measurement delay problem, such as predictive filtering proposed by Sieberling et al. [14]. However, the prediction requires additional modeling and cannot predict disturbances. The final reference command for x2ref can be given by
(19)x2ref=x2f+Δx2.


### 3.2.3. Uncertainty Analysis

The model uncertainty stemmed from aerodynamic force in translational dynamic control loop is analyzed in this subsection. The change of aerodynamic coefficients is assumed to be primarily caused by the angle of attack (α) and sideslip angle (β) for simplification. The aerodynamic coefficient can be expressed as
(20)CD=C¯D0+CDαα,CC=C¯C0+CCββ,CL=C¯L0+CLαα.In (20), CD, CC, and CL represent the total drag coefficient, side force coefficient, and lift coefficient, respectively; C¯D0 represents the sum of drag coefficients without the part contributed by α; C¯C0 represents the sum of side force coefficients without the part contributed by β; C¯L0 represents the sum of lift coefficients without the part contributed by α; and CDα, CCβ, and CLα represent derivative coefficients about α and β. These coefficients vary with the flight condition and can be assumed as constants in the calculation of each control period. The methods to obtain these parameters involve mathematic estimation and CFD analysis, as described in Section 2. The errors between the estimate coefficients and accurate coefficients make the model uncertain and assume that the uncertainty is mainly caused by C¯D0, C¯C0, and C¯L0. As the work by Lu et al. [18] indicates, the kinematic roll angle (μ) can be calculated directly under the assumption of β=0, which largely simplifies the control equation. Additionally, coefficients C¯D0, C¯C0, and C¯L0 can be eliminated by adopting INDI, which evidently reduces model uncertainty. However, the UAV in this paper with direct force control owns great lateral maneuverability. Therefore, β cannot be approximated as zero, and μ is calculated through solving the control equation. As a result, coefficients C¯D0, C¯C0, and C¯L0 arise in the third column of g1a and cannot be eliminated. Nevertheless, no uncertainty exists in the thrust vector control matrix g1t.The application of INDI in the translational dynamics control loop has two primary advantages. First, INDI control law restrains the model uncertainty caused by aerodynamic force in control matrixg1a and isolates its impact on direct force incremental control. Second, INDI can decouple the control variables in g1x1,x2 and transfer the control equation into a linear form. Accordingly, the aerodynamic force and vectored thrust control allocation are simplified, the complex numerical solution of nonlinear coupling control is avoided, and computational load of onboard control system is reduced.


## 3.2.1. INDI Control Law Design

In translational dynamic control loop, the vehicle’s flight path statusx1 is controlled by aerodynamic force and vectored thrust together to track the reference command x1ref obtained by the last control loop. Define the control variables as
(5)x2=αβμTxTyTzT,x2a=αβμT,x2t=TxTyTzT.In this loop, there are totally 6 control inputs. Among them,x2a represent the aerodynamic control variables, where μ is the bank angle, and α and β are the angle of attack and angle of sideslip, respectively. x2t represent the vectored thrusts, and Tx, Ty, and Tz are the components of the thrust in the body axis. The equation of translational dynamics is expressed as
(6)ẋ1=f1x1+g1x1,x2,(7)f1x1=1m0001mVcosγ000−1mVLkg00mg,(8)g1x1,x2=1m0001mVcosγ000−1mVLkbα,β,μx2t+LkaμFA.In (7) and (8), Lkg is the direction cosine matrix from earth axis to flight path axis, Lkbα,β,μ is the direction cosine matrix from body axis to flight path axis, and Lkaμ is the direction cosine matrix from wind axis to flight path axis. FA=−DαCβ−LαT is the aerodynamic force in which D represents drag, C represents aerodynamic side force, and L represents lift. FA and Lkbα,β,μ and Lkaμ are assumed as the functions of x2a. The multiplication of direction cosine matrixes causes severe coupling in control variables. Furthermore, in translational dynamic control loop, the error caused by the estimation of FA will bring uncertainties in the control system. Therefore, the INDI method is adopted to address the control problem. To rewrite (6) in incremental form, it is defined that x1n and x2n represent the values of x1 and x2 in the next time step, and their relationships are described below:
(9)x1n=x1+Δx1,x2n=x2+Δx2.Applying Taylor expansion toẋ1n at x1,x2 and with higher-order terms neglected, the result can be expressed as
(10)ẋ1n=ẋ1+∂f1x1∂x1x1n−x1+∂g1x1,x2∂x1x1n−x1+∂g1x1,x2∂x2x2n−x2.In (10), the second and third terms partial to x1 are assumed much smaller than the forth term, partial to x2. This commonly arises from the principle of time scale separation [20]. For simplification, it is approximated that
(11)∂g1x1,x2∂x2x2n−x2=g1x1,x2+Δx2−g1x1,x2Δx2x2n−x2=g1x1,x2+Δx2−g1x1,x2.Replaceẋ1n with ideal derivatives of flight path vector ẋ1des, which is calculated by the linear controller. The incremental control equation can be denoted as
(12)ẋ1des−ẋ1=g1x1,x2+Δx2−g1x1,x2.To simplify the computation process, it is assumed that(13)Δx2iΔx2j=0,where x2i and x2j represent the elements in vector x2. Accordingly, the equation can be rewritten into the affine-in-control form
(14)ẋ1des−ẋ1=g1a00g1tΔx2aΔx2t,where g1a and g1t are the 3×3 matrixes, representing the aerodynamic control matrix and thrust vector control matrix. The detail information of g1a and g1t is shown in Appendix. Compared with a traditional aircraft in incremental approach, the thrust vector control matrix g1t is additional. The elements in g1t are polynomials of control variables in x2a. Thus, the control system is transferred into a linear and time variant system, and the increment of control variables is decoupled. Then the ideal incremental virtual command ∆x2 can be calculated by
(15)Δx2aΔx2t=g1a00g1t‡ẋ1des−ẋ1.In (15), the superscript “‡” represents the dynamic weight pseudoinverse method, which will be introduced in the next section. The reference command for the next control loop is obtained by
(16)x2ref=x2ax2t+Δx2aΔx2t.


## 3.2.2. Flight Path Vector Derivative Acquisition

In (12), (14), and (15), ẋ1 can be derived by V, χ, and γ, which are measured by onboard sensors. With the ẋ1des−ẋ1 introduced, the term f1x1 is cancelled, which is the reason why INDI is referred to a sensor-based approach. This approach transfers the dependence on the model accuracy into the dependence on the sensor accuracy. In most cases, the signals were measured by sensor contain noise, and the differentiation of noisy signal amplifies the noise. To make results accurate, a filter can be adopted to abate the noise in sensor data. In this paper, a second-order filter is employed. As stated in literatures [16, 18, 21], the washout filter can be expressed in Laplace domain as
(17)Hs=ωn2s2+2ζωns+ωn2.However, the filter leads to a delay which should be compensated. In the Taylor expansion shown in (10), ẋ1, x1, and x2 should be from the same moment. In this regard, a second-order filter is also applied for x2 to counteract the impact caused by time delay in x1 and ẋ1, and (12) can be rewritten in a time-synchronized form
(18)ẋ1des−ẋ1f=g1x1f,x2f+Δx2−g1x1f,x2f,where subscript “f” represents the filtered variable. Alternatively, other methods can also deal with the measurement delay problem, such as predictive filtering proposed by Sieberling et al. [14]. However, the prediction requires additional modeling and cannot predict disturbances. The final reference command for x2ref can be given by
(19)x2ref=x2f+Δx2.


## 3.2.3. Uncertainty Analysis

The model uncertainty stemmed from aerodynamic force in translational dynamic control loop is analyzed in this subsection. The change of aerodynamic coefficients is assumed to be primarily caused by the angle of attack (α) and sideslip angle (β) for simplification. The aerodynamic coefficient can be expressed as
(20)CD=C¯D0+CDαα,CC=C¯C0+CCββ,CL=C¯L0+CLαα.In (20), CD, CC, and CL represent the total drag coefficient, side force coefficient, and lift coefficient, respectively; C¯D0 represents the sum of drag coefficients without the part contributed by α; C¯C0 represents the sum of side force coefficients without the part contributed by β; C¯L0 represents the sum of lift coefficients without the part contributed by α; and CDα, CCβ, and CLα represent derivative coefficients about α and β. These coefficients vary with the flight condition and can be assumed as constants in the calculation of each control period. The methods to obtain these parameters involve mathematic estimation and CFD analysis, as described in Section 2. The errors between the estimate coefficients and accurate coefficients make the model uncertain and assume that the uncertainty is mainly caused by C¯D0, C¯C0, and C¯L0. As the work by Lu et al. [18] indicates, the kinematic roll angle (μ) can be calculated directly under the assumption of β=0, which largely simplifies the control equation. Additionally, coefficients C¯D0, C¯C0, and C¯L0 can be eliminated by adopting INDI, which evidently reduces model uncertainty. However, the UAV in this paper with direct force control owns great lateral maneuverability. Therefore, β cannot be approximated as zero, and μ is calculated through solving the control equation. As a result, coefficients C¯D0, C¯C0, and C¯L0 arise in the third column of g1a and cannot be eliminated. Nevertheless, no uncertainty exists in the thrust vector control matrix g1t.The application of INDI in the translational dynamics control loop has two primary advantages. First, INDI control law restrains the model uncertainty caused by aerodynamic force in control matrixg1a and isolates its impact on direct force incremental control. Second, INDI can decouple the control variables in g1x1,x2 and transfer the control equation into a linear form. Accordingly, the aerodynamic force and vectored thrust control allocation are simplified, the complex numerical solution of nonlinear coupling control is avoided, and computational load of onboard control system is reduced.


## 3.3. Rotational Kinematic Control Loop

The objective in this control loop is to trackx2aref given in the translational dynamics control loop. The angular rate vector of the vehicle is defined as
(21)x3=pqrT.In (21), p, q, and r represent the roll, pitch, and yaw rates in body axis, respectively. No model uncertainty exists in this control loop, and the control law based on the standard NDI approach is established below:
(22)ẋ2a=f2x2a+g2x2ax3,f2x2a=1cosβ−cosγ−sinμcosγ−sinμcosβcosμcosγcosβcosμsinβsinγcosβ+sinμcosγsinβγ̇χ̇,g2x2a=1cosβ−cosαsinβcosβ−sinαsinβsinαcosβ0−cosαcosβcosα0sinα.It is noteworthy thatγ̇, χ̇, and μ are not directly measured onboard in this control system. They are calculated by (23) derived from (6).
(23)γ̇=cosαsinβsinμ+sinαcosμAx−cosβsinμAy+sinαsinβsinμ−cosαcosμAzV,χ̇=−cosαsinβcosμ+sinαsinμAx+cosβcosμAy+−sinαsinβcosμ−cosαsinμAzVcosγ,where Ax, Ay, and Az are directly measured by the accelerometers of the aircraft. The reason why γ̇ and χ̇ are not acquired by the same method adopted in translational dynamic control loop is that the filter will delay the measurement and cause errors from accurate results. The reference angular rates can be obtained by
(24)x3ref=g2x2a−1ẋ2ades−f2x2a,where ẋ2ades is the desired command and designed by linear controller, which is similar to the one adopted to design ẋ1des.


## 3.4. Rotational Dynamic Control Loop

x4 is defined as the control moments acting on the vehicle, and it can be denoted as
(25)x4=x4s+x4t=lcmcncT,where x4s denote the moments generated by aerodynamic control surface; x4t indicate the moments generated by vectored thrust; and lc, mc, and nc represent the rolling, pitching, and yawing control moments, respectively. The dynamics of the angular rates of the vehicle can be expressed into the following affine-in-control form
(26)ẋ3=f3x3+J−1x4,f3x3=J−1Ma−x3×Jx3,where J denotes the inertia matrix and Ma refers to the aerodynamic moments generated by derivatives unrelated to control surface deflections.
(27)Ma=q¯SC¯l0C¯m0C¯n0.Like aerodynamic force, the aerodynamic coefficients can also be denoted as(28)Cl=C¯l0+Clδaδa+Clδrδr,Cm=C¯m0+Cmδeδe,Cn=C¯n0+Cnδaδa+Cnδrδr.In (27) and (28), Cl, Cm, and Cn represent the total 3-axis aerodynamic moment coefficients, respectively; C¯l0 represents the sum of aerodynamic rolling moment coefficients without the part contributed by δa and δr; C¯m0 represents the sum of aerodynamic pitching moment coefficients without the part contributed by δe; C¯n0 represents the sum of aerodynamic yawing moment without the part contributed by δa and δr; and Clδa, Clδr, Cmδe, Cnδa, and Cnδr represent derivative aerodynamic moment coefficients about δa, δe, and δr.The aerodynamic control moments are denoted as(29)x4s=q¯Sb1CASδAS,CAS=Clδa0Clδr0Cmδe0Cnδa0Cnδr,b1=diagb,c¯,b,where q¯ denotes the dynamic pressure, S denotes the reference wing area, b denotes the wing span, and c¯ denotes the mean aerodynamic chord. δAS=δaδeδrT are control surface deflections. The control moment coefficients CAS are assumed to be accurate in aerodynamic estimation, and the uncertainties are considered to be caused principally by C¯l0, C¯m0, and C¯n0. Rewriting (26) into the incremental form, as presented in (30), the term f3x3 can be cleared up.
(30)Δx4=Jẋ3des−ẋ3f.In (30), ẋ3f represents the derivative of x3, which is the filtered signal collected by the sensor. In the absence of f3x3, the uncertainty caused by Ma is eliminated and the nonlinear cross couplings of the angular rate term is also cancelled. The reference control moments can be calculated as
(31)x4ref=x4f+Δx4.After the calculation of four control loops,∆x2t and ∆x4 will be output as virtual commands to the second-layer control allocation. The control allocation of the engine thrust, vectoring nozzles, and aerodynamic control surface deflections will be introduced in the next section. The block diagram of the control system is illustrated in Figure 2.Figure 2
Conceptual block diagram of the control system.


## 4. Two-Layer Cascaded Optimal Control Allocation

The allocation method of redundant control variables (elements inx2 and x4) in translational dynamics and rotational dynamics control loops is introduced in this section. Given this, the allocation result of x2 will affect x4, and the allocation results of power system actuators TL, TR, TF δR, δL,andδF are determined by x2t and x4t synthetically. A two-layer cascade optimal control allocation method is designed to address UAV’s control allocation based on INDI control law.


### 4.1. First-Layer Trajectory Incremental Control Allocation

The first-layer control allocation is conducted to allocate the increment of flight attitudeΔx2a and vectored thrust Δx2t in trajectory control.


#### 4.1.1. Incremental Pseudoinverse Method

In translational dynamics control loop, the control equation is linearized by the INDI method. Accordingly, the dynamic weight pseudoinverse method can be used in control allocation. Based on (14), it is defined as
(32)B=g1a00g1t,ν=ẋ1des−ẋ1f,where ẋ1f represents the derivative of x1, which is the filtered signal collected by the sensor. As the filter will delay the time in ẋ1f, to keep all signals synchronized in control equations, the signal of x2 is also required to be filtered to counteract the error generated by the delay of ẋ1f. Regardless of the control variable constrains, the optimal control allocation can be denoted as
(33)minΔx2W1Δx2+x2f2+W2Δx22(34)s.t.BΔx2=ν,where Δx2 denotes the optimization variable and x2f and ν are considered constants, they are all 3×1 vectors, and W1 and W2 indicate the diagonal matrixes of weight, which are generated by dynamic weight strategy, and will be discussed later. The first term of objective function, as shown in (33), is adopted to control the scale of control variables in x2, while the second term is employed to limit the change rates of the control variables. Equation (33) can be expanded as
(35)W1Δx2+x2f2+W2Δx22=Δx2+x2fTW12Δx2+x2f+xTW22Δx2=Δx2TW12Δx2+Δx2TW12x2f+x2fTW12Δx2+x2fTW12x2f+Δx2TW22Δx2.In (35), Δx2TW12x2f=x2fTW12Δx2 and x2fTW12x2f are constants which can be calculated directly. Accordingly, it is deduced that (33) has the same minimizing argument as
(36)Δx2TW12Δx2+2Δx2TW12x2f+Δx2TW22Δx2.Define(37)W=W12+W22,x20=W−2W12x2f.Substitute (37) into (36) and incorporate with constant x20TW2x20.
(38)Δx2TW2Δx2+2Δx2TW2x20+x20TW2x20=Δx2+x20TW2Δx2+x20=WΔx2+x202.Based on the foregoing deduction, the optimal problem presented in (33) is equivalence
(39)minΔx2WΔx2+x202s.t.BΔx2=ν.On that basis, the minimum norm solution of the control allocation problem is obtained as(40)Δx2=Gν+GB−Ix20,(41)G=W−1BW−1†.In (41), superscript “†” represents the pseudoinverse operation for the matrix. In order to make the first-layer control allocation an equality constrained optimization problem, the range of Δx2 is not limited. Then an analytical solution is obtained, and numerical iterations are avoided, which is suitable for onboard use [22]. With the allocation results of Δx2a, the x2aref can be calculated by (19). After adjusting x2aref into the right quadrant, as shown in (42), x2aref is transmitted into the rotational kinematics control loop.
(42)x2aref∈−0.5π0.5π.TheΔx2t is considered the virtual command to transmit into the second-layer control allocation, in which Δx2t and Δx4t are used to solve the values of power system actuators TL, TR, TF δR, δL,andδF.


#### 4.1.2. Dynamic Weight Strategy

In the course of control allocation, the weight is adopted to describe the differences of control variables’ significance. On the basis of the significance of control variables changing with flight states and mission requirements, a dynamic weight strategy is proposed to allocate control variables optimally and properly.The traditional determination of control variable’s weight largely depends on human experience. In this regard, there is no absolute criteria for weight determination. It is commonly difficult to judge a control variable’s significance globally, while the significance between every two control variables can be easily compared. In dynamic weight strategy, the analytic hierarchy process (AHP) is adopted to synthesize comparison results of every two control variables and calculate each control variable’s weight by AHP-judgment matrix. The weight matrixes in the objective function can be denoted as(43)W1=diagWα,Wβ,Wμ,WTx,WTy,WTz,W2=diagWΔα,WΔβ,WΔμ,WΔTx,WΔTy,WΔTz.In (43), each element in W1 and W2 represents the weight of the control variable corresponding to its subscript. According to different physical properties, control variables are divided into different sets. The definitions and subordinations, respectively, of the sets are expressed in
(44)Sx2=α,β,μ,Tx,Ty,Tz,Sx2a=α,β,μ,Sx2t=Tx,Ty,Tz,SΔx2=Δα,Δβ,Δμ,ΔTx,ΔTy,ΔTz,SΔx2a=Δα,Δβ,Δμ,SΔx2t=ΔTx,ΔTy,ΔTz,(45)Sx2a⊆Sx2,Sx2t⊆Sx2,SΔx2a⊆SΔx2,SΔx2t⊆SΔx2.Based on the subordinations, the sets and control variables are classified into three hierarchies, known as the weight structure, as presented in Figure3.Figure 3
Illustration of weight structure.The weights in different hierarchies are determined, respectively, and they all obey the rules listed below:(1)
Weights in the same set are required to be normalized(46)∑i=12wi=1,∑j=12wij=1i=1,2,∑k=13wijk=1i,j=1,2
(2)
Weights in the first and second hierarchies satisfy(47)w1,w2,w11,w12,w21,w22∈0:0.01:1
(3)
In the third hierarchy,wijk is acquired by the APH-judgment matrix. The judgment matrix is established for each set, with elements representing the priorities of every two control variables in the same set. The processes of APH method are introduced in [23, 24]. According to the weights in each hierarchy, the final weights in the objective function can be calculated by
(48)W1k,k=w1×w11×w11kk=1,2,3,W1k+3,k+3=w1×w12×w12kk=1,2,3,W2k,k=w2×w21×w21kk=1,2,3,W2k+3,k+3=w2×w22×w22kk=1,2,3,where Wik,k represents the element in weight matrix Wi row k column k. In practical flight, the controller is required to generate the initial value of the weight structure according to the mission requirement and dynamically adjust the value of the weight structure according to the conditions of actuator saturation. In dynamic weight strategy, the weight generator and weight regulator are designed to realize the foregoing functions.The weight generator is designed to generate the initial value of the weight structure, in every control period before optimization. Its working principle is introduced below:Step 1.
Establish weight structures which consist ofwi, wij, and wijk.Step 2.
Artificially design and test several weight structures according to some typical flight states and mission requirements. On this basis, index weight structures with their relative flight states and mission requirements save weight structures into repository.Step 3.
Acquire weight structure corresponding to the current flight state and mission requirement through traversing the repository.A simplified repository is established in this paper. For flight states, merely the impact of velocity is factored in, and velocities 5 m/s, 10 m/s, 15 m/s, 20 m/s, and 25 m/s count as typical flight states. In mission requirement, only the impact of direct force control level (DFClevel∈0,1) is factored in, and DFClevel equals 0, 0.5, and 1 count as typical mission requirements. Larger DFClevel represents less attitude changes and more direct force control, while smaller DFClevel which represents less vectored thrust usage and more attitude maneuver during trajectory tracking. The repository of 3 × 5 weight structures employed in this paper is acquired by artificial design and computer-based simulation. In control allocation, the weight structure for the current demand is interpolated by velocity and DFClevel.The weight regulator is designed to adjust the weight when control variables exceed their limitation and ensure the rationality of allocation results. The working process of weight regulator is shown below.Step 1.
Extract control variables saturated in the last control allocation according to the feedback information. The saturated control variable can be single or multiple.Step 2.
Update saturation counter. Every control variable has a corresponding saturation counternijk∈Ni,j=1,2,k=1,2,3, which is adopted to record saturation times. If the control variable is not saturated in the last allocation, its saturation counter will be established at zero.Step 3.
Adjust the weight structure in accordance with the saturation counter results. For different hierarchies of weight structure, the adjustment strategies are different, as illustrated below:Strategy 1.
Fornijk∈1,150 with the interval of every 5 count times, increase the significance of the saturated control variables relative to other control variables in the identical minimum set. On that basis, update the APH-judgment matrix and recalculate wijk in set. If more than one control variable gets saturated, the comparison between the saturated control variables in significance remains unchanged, and their significance is improved compared with those unsaturated.Strategy 2.
Fornijk∈151,300, the weight (wij) in the middle hierarchy is required to be adjusted as Strategy 1 is exercised. Take the angle of attack α as an example. If α continues to be saturated with the interval of every 5 count times, increase the weight (w11) of set Sx2a. In (47), wij is divided into 101 levels, and wij is improved to a new level for each weight adjustment. Assume w11=0.6andw12=0.4 before the adjustment, and with the improvement of w11, the values turn out to be w11=0.61andw12=0.39.Strategy 3.
Fornijk∈301,+∞ with an interval of every 10 count times, the weight (wi) of set (Sx2 or S∆x2) is required to be adjusted, as Strategy 1 and Strategy 2 are exercised. If the weight reaches the upper limit with wi,wij=1, the weight stops increasing; if the weight reaches the lower limit with wi,wij=0, the decrement stops.The range of saturation counternijk and the adjustment of interval in each strategy are determined by control frequency in the weight regulator. For example, for 100 Hz control frequency, the control is allocated for every 0.01 second. Following the foregoing strategy, if a control variable is saturated for more than 3 seconds, the controller will adjust the weight for every hierarchy in the weight structure to get optimal control allocation to satisfy the current flight requirement. For the excessively large range of nijk, the weight adjustment turns out to be over-sluggish. The control variables will be saturated continuously, and the vehicle will deviate from target trajectory seriously. For the excessively small range of nijk, the weight adjustment will be excessively sensitive. The fast-changing weight in objective function will cause vehicle oscillation and control divergence. The working process of the first-layer control allocation and dynamic weight strategy is illustrated in Figure 4.Figure 4
Working process of first-layer control allocation and dynamic weight strategy.


### 4.2. Second-Layer Actuator Control Allocation

In the second-layer control allocation, the allocation of aerodynamic control surfaces (δa, δe, δr) and power system actuators (TL, TR, TF, δL, δR, δF) are conducted according to the information of vectored thrust and control moments calculated in translational dynamics control loop and rotational dynamics control loops, respectively.The dynamic model of the engine system can be denoted as(49)x2t=TxTyTz=TLcosδL+TRcosδRTFsinδFTLsinδL+TRsinδR−TFcosδF,x4t=lTmTnT=−TLsinδLdyt+TRsinδRdytTLsinδL+TLsinδLdxt+TFcosδFdxfTLcosδL−TRcosδRdyt+TFsinδFdxf.Define(50)TE=TLTRTFT,δTVN=δLδRδFT,δAS=δaδeδrT.In (49) and (50), TF denotes lift fan thrust; TL and TR denote the left and right engine thrusts, respectively; δF denotes the lift fan rudder deflection angle; δL and δR denote the left and right engine thrust vector nozzle deflections, respectively; dxf denotes the distance from the lift fan center to the vehicle mass center in body x-axis; dxt and dyt denote the distance from the engine nozzle to the vehicle mass center in body x-axis and y-axis, respectively; and δa, δe, and δr represent the alerion, elevator, and rudder, respectively. To keep time synchronized, the output of the actuator model will go through a second-order filter to get TE.f,δTVN.f,δAS.f, where subscript “f” denotes the filtered control output. In the second-layer control allocation, the daisy chaining method is adopted to allocate aerodynamic control moment and thrust vector moment. On that basis, the control output of the power system is calculated by x2t and x4t. The second-layer control allocation can be conducted in the incremental linear form or normal nonlinear form, and the two methods are introduced below.


#### 4.2.1. Incremental Linear Allocation Method

The incremental linear allocation method allocates∆TE ∆δTVN and ∆δAS according to ∆x4 and ∆x2t. The incremental daisy chaining is adopted to allocate ∆x4s and ∆x4t based on (51). In this method, the priority of ∆x4s outstrips ∆x4t, and thus, the thrust vector moment is adopted only when the aerodynamic control moment gets saturated.
(51)Δx4=Δx4s+Δx4t.The increment of aerodynamic control surface is calculated by(52)ΔδAS=Δx4sq¯Sb1CAS.The reference value of aerodynamic control surface can be expressed as(53)δASref=δAS.f+ΔδAS.On that basis, introduce (49) into the incremental form, as shown in (54), in which the incremental multiplication is considered a high-order small quantity and can be ignored.
(54)Δx2tΔx4t=GEΔTEΔδTVN,(55)GE=cosδLcosδR0−TLsinδL−TRsinδR000sinδF00TFcosδFsinδLsinδR−cosδFTLcosδLTRcosδRTFsinδF−sinδLdytsinδRdyt0−TLcosδLdytTRcosδRdyt0sinδLdxtsinδRdxtcosδFdxfTLcosδLdxtTRcosδRdxt−TFsinδFdxfcosδLdyt−cosδRdytsinδFdxf−TLsinδLdytTRsinδRdytTFcosδFdxf.The thrust vector engine model is transformed from a nonlinear system into a linear system, and∆TE and ∆δTVN can be calculated directly by matrix inversion through the incremental linear allocation, as shown in (56). This method can evidently expedite allocation without iterative solution of nonlinear equations. However, like (14), when UAV is in high maneuver flight, the increment of control variables will be large, and the neglection of incremental multiplication will cause approximation errors. In this regard, this method is appropriate for minor or moderate maneuverable flights.
(56)ΔTEΔδTVN=GE−1Δx2tΔx4t.The reference thrust vector engine actuators are calculated by (57).
(57)TEref=TE.f+ΔTE,δTVNref=δTVN.f+ΔδTVN.The range ofTL, TR, TF and δL, δR, δF is limited on the basis of power system, as shown in (58). Eventually, the control signals are output to the aerodynamic control surfaces and power system.
(58)Ti=Timax,ifTi≥Timax,Ti=Timin,ifTi≤Timin,δi=minδitn−1+δ̇imaxΔt,δimax,ifδitn≥minδitn−1+δ̇imaxΔt,δimax,δi=maxδitn−1−δ̇imaxΔt,δimin,ifδitn≤maxδitn−1−δ̇imaxΔt,δimin.The flow chart of the incremental linear allocation is presented in Figure5.Figure 5
Second-layer incremental linear allocation method.


#### 4.2.2. Normal Nonlinear Allocation Method

In nonlinear allocation, other than calculating∆TE and ∆δTVN, the outputs of thrust vector engine system TEandδTVN are calculated directly through solving nonlinear equations. This method is suitable for all maneuver flights and will not cause approximation error. The total time consumption remains rational and applicable for onboard use though the calculation time is longer than the incremental linear allocation method. The flow chart of the allocation process is exhibited in Figure 6.Figure 6
Second-layer normal nonlinear allocation method.


## 4.1. First-Layer Trajectory Incremental Control Allocation

The first-layer control allocation is conducted to allocate the increment of flight attitudeΔx2a and vectored thrust Δx2t in trajectory control.


### 4.1.1. Incremental Pseudoinverse Method

In translational dynamics control loop, the control equation is linearized by the INDI method. Accordingly, the dynamic weight pseudoinverse method can be used in control allocation. Based on (14), it is defined as
(32)B=g1a00g1t,ν=ẋ1des−ẋ1f,where ẋ1f represents the derivative of x1, which is the filtered signal collected by the sensor. As the filter will delay the time in ẋ1f, to keep all signals synchronized in control equations, the signal of x2 is also required to be filtered to counteract the error generated by the delay of ẋ1f. Regardless of the control variable constrains, the optimal control allocation can be denoted as
(33)minΔx2W1Δx2+x2f2+W2Δx22(34)s.t.BΔx2=ν,where Δx2 denotes the optimization variable and x2f and ν are considered constants, they are all 3×1 vectors, and W1 and W2 indicate the diagonal matrixes of weight, which are generated by dynamic weight strategy, and will be discussed later. The first term of objective function, as shown in (33), is adopted to control the scale of control variables in x2, while the second term is employed to limit the change rates of the control variables. Equation (33) can be expanded as
(35)W1Δx2+x2f2+W2Δx22=Δx2+x2fTW12Δx2+x2f+xTW22Δx2=Δx2TW12Δx2+Δx2TW12x2f+x2fTW12Δx2+x2fTW12x2f+Δx2TW22Δx2.In (35), Δx2TW12x2f=x2fTW12Δx2 and x2fTW12x2f are constants which can be calculated directly. Accordingly, it is deduced that (33) has the same minimizing argument as
(36)Δx2TW12Δx2+2Δx2TW12x2f+Δx2TW22Δx2.Define(37)W=W12+W22,x20=W−2W12x2f.Substitute (37) into (36) and incorporate with constant x20TW2x20.
(38)Δx2TW2Δx2+2Δx2TW2x20+x20TW2x20=Δx2+x20TW2Δx2+x20=WΔx2+x202.Based on the foregoing deduction, the optimal problem presented in (33) is equivalence
(39)minΔx2WΔx2+x202s.t.BΔx2=ν.On that basis, the minimum norm solution of the control allocation problem is obtained as(40)Δx2=Gν+GB−Ix20,(41)G=W−1BW−1†.In (41), superscript “†” represents the pseudoinverse operation for the matrix. In order to make the first-layer control allocation an equality constrained optimization problem, the range of Δx2 is not limited. Then an analytical solution is obtained, and numerical iterations are avoided, which is suitable for onboard use [22]. With the allocation results of Δx2a, the x2aref can be calculated by (19). After adjusting x2aref into the right quadrant, as shown in (42), x2aref is transmitted into the rotational kinematics control loop.
(42)x2aref∈−0.5π0.5π.TheΔx2t is considered the virtual command to transmit into the second-layer control allocation, in which Δx2t and Δx4t are used to solve the values of power system actuators TL, TR, TF δR, δL,andδF.


### 4.1.2. Dynamic Weight Strategy

In the course of control allocation, the weight is adopted to describe the differences of control variables’ significance. On the basis of the significance of control variables changing with flight states and mission requirements, a dynamic weight strategy is proposed to allocate control variables optimally and properly.The traditional determination of control variable’s weight largely depends on human experience. In this regard, there is no absolute criteria for weight determination. It is commonly difficult to judge a control variable’s significance globally, while the significance between every two control variables can be easily compared. In dynamic weight strategy, the analytic hierarchy process (AHP) is adopted to synthesize comparison results of every two control variables and calculate each control variable’s weight by AHP-judgment matrix. The weight matrixes in the objective function can be denoted as(43)W1=diagWα,Wβ,Wμ,WTx,WTy,WTz,W2=diagWΔα,WΔβ,WΔμ,WΔTx,WΔTy,WΔTz.In (43), each element in W1 and W2 represents the weight of the control variable corresponding to its subscript. According to different physical properties, control variables are divided into different sets. The definitions and subordinations, respectively, of the sets are expressed in
(44)Sx2=α,β,μ,Tx,Ty,Tz,Sx2a=α,β,μ,Sx2t=Tx,Ty,Tz,SΔx2=Δα,Δβ,Δμ,ΔTx,ΔTy,ΔTz,SΔx2a=Δα,Δβ,Δμ,SΔx2t=ΔTx,ΔTy,ΔTz,(45)Sx2a⊆Sx2,Sx2t⊆Sx2,SΔx2a⊆SΔx2,SΔx2t⊆SΔx2.Based on the subordinations, the sets and control variables are classified into three hierarchies, known as the weight structure, as presented in Figure3.Figure 3
Illustration of weight structure.The weights in different hierarchies are determined, respectively, and they all obey the rules listed below:(1)
Weights in the same set are required to be normalized(46)∑i=12wi=1,∑j=12wij=1i=1,2,∑k=13wijk=1i,j=1,2
(2)
Weights in the first and second hierarchies satisfy(47)w1,w2,w11,w12,w21,w22∈0:0.01:1
(3)
In the third hierarchy,wijk is acquired by the APH-judgment matrix. The judgment matrix is established for each set, with elements representing the priorities of every two control variables in the same set. The processes of APH method are introduced in [23, 24]. According to the weights in each hierarchy, the final weights in the objective function can be calculated by
(48)W1k,k=w1×w11×w11kk=1,2,3,W1k+3,k+3=w1×w12×w12kk=1,2,3,W2k,k=w2×w21×w21kk=1,2,3,W2k+3,k+3=w2×w22×w22kk=1,2,3,where Wik,k represents the element in weight matrix Wi row k column k. In practical flight, the controller is required to generate the initial value of the weight structure according to the mission requirement and dynamically adjust the value of the weight structure according to the conditions of actuator saturation. In dynamic weight strategy, the weight generator and weight regulator are designed to realize the foregoing functions.The weight generator is designed to generate the initial value of the weight structure, in every control period before optimization. Its working principle is introduced below:Step 1.
Establish weight structures which consist ofwi, wij, and wijk.Step 2.
Artificially design and test several weight structures according to some typical flight states and mission requirements. On this basis, index weight structures with their relative flight states and mission requirements save weight structures into repository.Step 3.
Acquire weight structure corresponding to the current flight state and mission requirement through traversing the repository.A simplified repository is established in this paper. For flight states, merely the impact of velocity is factored in, and velocities 5 m/s, 10 m/s, 15 m/s, 20 m/s, and 25 m/s count as typical flight states. In mission requirement, only the impact of direct force control level (DFClevel∈0,1) is factored in, and DFClevel equals 0, 0.5, and 1 count as typical mission requirements. Larger DFClevel represents less attitude changes and more direct force control, while smaller DFClevel which represents less vectored thrust usage and more attitude maneuver during trajectory tracking. The repository of 3 × 5 weight structures employed in this paper is acquired by artificial design and computer-based simulation. In control allocation, the weight structure for the current demand is interpolated by velocity and DFClevel.The weight regulator is designed to adjust the weight when control variables exceed their limitation and ensure the rationality of allocation results. The working process of weight regulator is shown below.Step 1.
Extract control variables saturated in the last control allocation according to the feedback information. The saturated control variable can be single or multiple.Step 2.
Update saturation counter. Every control variable has a corresponding saturation counternijk∈Ni,j=1,2,k=1,2,3, which is adopted to record saturation times. If the control variable is not saturated in the last allocation, its saturation counter will be established at zero.Step 3.
Adjust the weight structure in accordance with the saturation counter results. For different hierarchies of weight structure, the adjustment strategies are different, as illustrated below:Strategy 1.
Fornijk∈1,150 with the interval of every 5 count times, increase the significance of the saturated control variables relative to other control variables in the identical minimum set. On that basis, update the APH-judgment matrix and recalculate wijk in set. If more than one control variable gets saturated, the comparison between the saturated control variables in significance remains unchanged, and their significance is improved compared with those unsaturated.Strategy 2.
Fornijk∈151,300, the weight (wij) in the middle hierarchy is required to be adjusted as Strategy 1 is exercised. Take the angle of attack α as an example. If α continues to be saturated with the interval of every 5 count times, increase the weight (w11) of set Sx2a. In (47), wij is divided into 101 levels, and wij is improved to a new level for each weight adjustment. Assume w11=0.6andw12=0.4 before the adjustment, and with the improvement of w11, the values turn out to be w11=0.61andw12=0.39.Strategy 3.
Fornijk∈301,+∞ with an interval of every 10 count times, the weight (wi) of set (Sx2 or S∆x2) is required to be adjusted, as Strategy 1 and Strategy 2 are exercised. If the weight reaches the upper limit with wi,wij=1, the weight stops increasing; if the weight reaches the lower limit with wi,wij=0, the decrement stops.The range of saturation counternijk and the adjustment of interval in each strategy are determined by control frequency in the weight regulator. For example, for 100 Hz control frequency, the control is allocated for every 0.01 second. Following the foregoing strategy, if a control variable is saturated for more than 3 seconds, the controller will adjust the weight for every hierarchy in the weight structure to get optimal control allocation to satisfy the current flight requirement. For the excessively large range of nijk, the weight adjustment turns out to be over-sluggish. The control variables will be saturated continuously, and the vehicle will deviate from target trajectory seriously. For the excessively small range of nijk, the weight adjustment will be excessively sensitive. The fast-changing weight in objective function will cause vehicle oscillation and control divergence. The working process of the first-layer control allocation and dynamic weight strategy is illustrated in Figure 4.Figure 4
Working process of first-layer control allocation and dynamic weight strategy.


## 4.1.1. Incremental Pseudoinverse Method

In translational dynamics control loop, the control equation is linearized by the INDI method. Accordingly, the dynamic weight pseudoinverse method can be used in control allocation. Based on (14), it is defined as
(32)B=g1a00g1t,ν=ẋ1des−ẋ1f,where ẋ1f represents the derivative of x1, which is the filtered signal collected by the sensor. As the filter will delay the time in ẋ1f, to keep all signals synchronized in control equations, the signal of x2 is also required to be filtered to counteract the error generated by the delay of ẋ1f. Regardless of the control variable constrains, the optimal control allocation can be denoted as
(33)minΔx2W1Δx2+x2f2+W2Δx22(34)s.t.BΔx2=ν,where Δx2 denotes the optimization variable and x2f and ν are considered constants, they are all 3×1 vectors, and W1 and W2 indicate the diagonal matrixes of weight, which are generated by dynamic weight strategy, and will be discussed later. The first term of objective function, as shown in (33), is adopted to control the scale of control variables in x2, while the second term is employed to limit the change rates of the control variables. Equation (33) can be expanded as
(35)W1Δx2+x2f2+W2Δx22=Δx2+x2fTW12Δx2+x2f+xTW22Δx2=Δx2TW12Δx2+Δx2TW12x2f+x2fTW12Δx2+x2fTW12x2f+Δx2TW22Δx2.In (35), Δx2TW12x2f=x2fTW12Δx2 and x2fTW12x2f are constants which can be calculated directly. Accordingly, it is deduced that (33) has the same minimizing argument as
(36)Δx2TW12Δx2+2Δx2TW12x2f+Δx2TW22Δx2.Define(37)W=W12+W22,x20=W−2W12x2f.Substitute (37) into (36) and incorporate with constant x20TW2x20.
(38)Δx2TW2Δx2+2Δx2TW2x20+x20TW2x20=Δx2+x20TW2Δx2+x20=WΔx2+x202.Based on the foregoing deduction, the optimal problem presented in (33) is equivalence
(39)minΔx2WΔx2+x202s.t.BΔx2=ν.On that basis, the minimum norm solution of the control allocation problem is obtained as(40)Δx2=Gν+GB−Ix20,(41)G=W−1BW−1†.In (41), superscript “†” represents the pseudoinverse operation for the matrix. In order to make the first-layer control allocation an equality constrained optimization problem, the range of Δx2 is not limited. Then an analytical solution is obtained, and numerical iterations are avoided, which is suitable for onboard use [22]. With the allocation results of Δx2a, the x2aref can be calculated by (19). After adjusting x2aref into the right quadrant, as shown in (42), x2aref is transmitted into the rotational kinematics control loop.
(42)x2aref∈−0.5π0.5π.TheΔx2t is considered the virtual command to transmit into the second-layer control allocation, in which Δx2t and Δx4t are used to solve the values of power system actuators TL, TR, TF δR, δL,andδF.


## 4.1.2. Dynamic Weight Strategy

In the course of control allocation, the weight is adopted to describe the differences of control variables’ significance. On the basis of the significance of control variables changing with flight states and mission requirements, a dynamic weight strategy is proposed to allocate control variables optimally and properly.The traditional determination of control variable’s weight largely depends on human experience. In this regard, there is no absolute criteria for weight determination. It is commonly difficult to judge a control variable’s significance globally, while the significance between every two control variables can be easily compared. In dynamic weight strategy, the analytic hierarchy process (AHP) is adopted to synthesize comparison results of every two control variables and calculate each control variable’s weight by AHP-judgment matrix. The weight matrixes in the objective function can be denoted as(43)W1=diagWα,Wβ,Wμ,WTx,WTy,WTz,W2=diagWΔα,WΔβ,WΔμ,WΔTx,WΔTy,WΔTz.In (43), each element in W1 and W2 represents the weight of the control variable corresponding to its subscript. According to different physical properties, control variables are divided into different sets. The definitions and subordinations, respectively, of the sets are expressed in
(44)Sx2=α,β,μ,Tx,Ty,Tz,Sx2a=α,β,μ,Sx2t=Tx,Ty,Tz,SΔx2=Δα,Δβ,Δμ,ΔTx,ΔTy,ΔTz,SΔx2a=Δα,Δβ,Δμ,SΔx2t=ΔTx,ΔTy,ΔTz,(45)Sx2a⊆Sx2,Sx2t⊆Sx2,SΔx2a⊆SΔx2,SΔx2t⊆SΔx2.Based on the subordinations, the sets and control variables are classified into three hierarchies, known as the weight structure, as presented in Figure3.Figure 3
Illustration of weight structure.The weights in different hierarchies are determined, respectively, and they all obey the rules listed below:(1)
Weights in the same set are required to be normalized(46)∑i=12wi=1,∑j=12wij=1i=1,2,∑k=13wijk=1i,j=1,2
(2)
Weights in the first and second hierarchies satisfy(47)w1,w2,w11,w12,w21,w22∈0:0.01:1
(3)
In the third hierarchy,wijk is acquired by the APH-judgment matrix. The judgment matrix is established for each set, with elements representing the priorities of every two control variables in the same set. The processes of APH method are introduced in [23, 24]. According to the weights in each hierarchy, the final weights in the objective function can be calculated by
(48)W1k,k=w1×w11×w11kk=1,2,3,W1k+3,k+3=w1×w12×w12kk=1,2,3,W2k,k=w2×w21×w21kk=1,2,3,W2k+3,k+3=w2×w22×w22kk=1,2,3,where Wik,k represents the element in weight matrix Wi row k column k. In practical flight, the controller is required to generate the initial value of the weight structure according to the mission requirement and dynamically adjust the value of the weight structure according to the conditions of actuator saturation. In dynamic weight strategy, the weight generator and weight regulator are designed to realize the foregoing functions.The weight generator is designed to generate the initial value of the weight structure, in every control period before optimization. Its working principle is introduced below:Step 1.
Establish weight structures which consist ofwi, wij, and wijk.Step 2.
Artificially design and test several weight structures according to some typical flight states and mission requirements. On this basis, index weight structures with their relative flight states and mission requirements save weight structures into repository.Step 3.
Acquire weight structure corresponding to the current flight state and mission requirement through traversing the repository.A simplified repository is established in this paper. For flight states, merely the impact of velocity is factored in, and velocities 5 m/s, 10 m/s, 15 m/s, 20 m/s, and 25 m/s count as typical flight states. In mission requirement, only the impact of direct force control level (DFClevel∈0,1) is factored in, and DFClevel equals 0, 0.5, and 1 count as typical mission requirements. Larger DFClevel represents less attitude changes and more direct force control, while smaller DFClevel which represents less vectored thrust usage and more attitude maneuver during trajectory tracking. The repository of 3 × 5 weight structures employed in this paper is acquired by artificial design and computer-based simulation. In control allocation, the weight structure for the current demand is interpolated by velocity and DFClevel.The weight regulator is designed to adjust the weight when control variables exceed their limitation and ensure the rationality of allocation results. The working process of weight regulator is shown below.Step 1.
Extract control variables saturated in the last control allocation according to the feedback information. The saturated control variable can be single or multiple.Step 2.
Update saturation counter. Every control variable has a corresponding saturation counternijk∈Ni,j=1,2,k=1,2,3, which is adopted to record saturation times. If the control variable is not saturated in the last allocation, its saturation counter will be established at zero.Step 3.
Adjust the weight structure in accordance with the saturation counter results. For different hierarchies of weight structure, the adjustment strategies are different, as illustrated below:Strategy 1.
Fornijk∈1,150 with the interval of every 5 count times, increase the significance of the saturated control variables relative to other control variables in the identical minimum set. On that basis, update the APH-judgment matrix and recalculate wijk in set. If more than one control variable gets saturated, the comparison between the saturated control variables in significance remains unchanged, and their significance is improved compared with those unsaturated.Strategy 2.
Fornijk∈151,300, the weight (wij) in the middle hierarchy is required to be adjusted as Strategy 1 is exercised. Take the angle of attack α as an example. If α continues to be saturated with the interval of every 5 count times, increase the weight (w11) of set Sx2a. In (47), wij is divided into 101 levels, and wij is improved to a new level for each weight adjustment. Assume w11=0.6andw12=0.4 before the adjustment, and with the improvement of w11, the values turn out to be w11=0.61andw12=0.39.Strategy 3.
Fornijk∈301,+∞ with an interval of every 10 count times, the weight (wi) of set (Sx2 or S∆x2) is required to be adjusted, as Strategy 1 and Strategy 2 are exercised. If the weight reaches the upper limit with wi,wij=1, the weight stops increasing; if the weight reaches the lower limit with wi,wij=0, the decrement stops.The range of saturation counternijk and the adjustment of interval in each strategy are determined by control frequency in the weight regulator. For example, for 100 Hz control frequency, the control is allocated for every 0.01 second. Following the foregoing strategy, if a control variable is saturated for more than 3 seconds, the controller will adjust the weight for every hierarchy in the weight structure to get optimal control allocation to satisfy the current flight requirement. For the excessively large range of nijk, the weight adjustment turns out to be over-sluggish. The control variables will be saturated continuously, and the vehicle will deviate from target trajectory seriously. For the excessively small range of nijk, the weight adjustment will be excessively sensitive. The fast-changing weight in objective function will cause vehicle oscillation and control divergence. The working process of the first-layer control allocation and dynamic weight strategy is illustrated in Figure 4.Figure 4
Working process of first-layer control allocation and dynamic weight strategy.


## 4.2. Second-Layer Actuator Control Allocation

In the second-layer control allocation, the allocation of aerodynamic control surfaces (δa, δe, δr) and power system actuators (TL, TR, TF, δL, δR, δF) are conducted according to the information of vectored thrust and control moments calculated in translational dynamics control loop and rotational dynamics control loops, respectively.The dynamic model of the engine system can be denoted as(49)x2t=TxTyTz=TLcosδL+TRcosδRTFsinδFTLsinδL+TRsinδR−TFcosδF,x4t=lTmTnT=−TLsinδLdyt+TRsinδRdytTLsinδL+TLsinδLdxt+TFcosδFdxfTLcosδL−TRcosδRdyt+TFsinδFdxf.Define(50)TE=TLTRTFT,δTVN=δLδRδFT,δAS=δaδeδrT.In (49) and (50), TF denotes lift fan thrust; TL and TR denote the left and right engine thrusts, respectively; δF denotes the lift fan rudder deflection angle; δL and δR denote the left and right engine thrust vector nozzle deflections, respectively; dxf denotes the distance from the lift fan center to the vehicle mass center in body x-axis; dxt and dyt denote the distance from the engine nozzle to the vehicle mass center in body x-axis and y-axis, respectively; and δa, δe, and δr represent the alerion, elevator, and rudder, respectively. To keep time synchronized, the output of the actuator model will go through a second-order filter to get TE.f,δTVN.f,δAS.f, where subscript “f” denotes the filtered control output. In the second-layer control allocation, the daisy chaining method is adopted to allocate aerodynamic control moment and thrust vector moment. On that basis, the control output of the power system is calculated by x2t and x4t. The second-layer control allocation can be conducted in the incremental linear form or normal nonlinear form, and the two methods are introduced below.


### 4.2.1. Incremental Linear Allocation Method

The incremental linear allocation method allocates∆TE ∆δTVN and ∆δAS according to ∆x4 and ∆x2t. The incremental daisy chaining is adopted to allocate ∆x4s and ∆x4t based on (51). In this method, the priority of ∆x4s outstrips ∆x4t, and thus, the thrust vector moment is adopted only when the aerodynamic control moment gets saturated.
(51)Δx4=Δx4s+Δx4t.The increment of aerodynamic control surface is calculated by(52)ΔδAS=Δx4sq¯Sb1CAS.The reference value of aerodynamic control surface can be expressed as(53)δASref=δAS.f+ΔδAS.On that basis, introduce (49) into the incremental form, as shown in (54), in which the incremental multiplication is considered a high-order small quantity and can be ignored.
(54)Δx2tΔx4t=GEΔTEΔδTVN,(55)GE=cosδLcosδR0−TLsinδL−TRsinδR000sinδF00TFcosδFsinδLsinδR−cosδFTLcosδLTRcosδRTFsinδF−sinδLdytsinδRdyt0−TLcosδLdytTRcosδRdyt0sinδLdxtsinδRdxtcosδFdxfTLcosδLdxtTRcosδRdxt−TFsinδFdxfcosδLdyt−cosδRdytsinδFdxf−TLsinδLdytTRsinδRdytTFcosδFdxf.The thrust vector engine model is transformed from a nonlinear system into a linear system, and∆TE and ∆δTVN can be calculated directly by matrix inversion through the incremental linear allocation, as shown in (56). This method can evidently expedite allocation without iterative solution of nonlinear equations. However, like (14), when UAV is in high maneuver flight, the increment of control variables will be large, and the neglection of incremental multiplication will cause approximation errors. In this regard, this method is appropriate for minor or moderate maneuverable flights.
(56)ΔTEΔδTVN=GE−1Δx2tΔx4t.The reference thrust vector engine actuators are calculated by (57).
(57)TEref=TE.f+ΔTE,δTVNref=δTVN.f+ΔδTVN.The range ofTL, TR, TF and δL, δR, δF is limited on the basis of power system, as shown in (58). Eventually, the control signals are output to the aerodynamic control surfaces and power system.
(58)Ti=Timax,ifTi≥Timax,Ti=Timin,ifTi≤Timin,δi=minδitn−1+δ̇imaxΔt,δimax,ifδitn≥minδitn−1+δ̇imaxΔt,δimax,δi=maxδitn−1−δ̇imaxΔt,δimin,ifδitn≤maxδitn−1−δ̇imaxΔt,δimin.The flow chart of the incremental linear allocation is presented in Figure5.Figure 5
Second-layer incremental linear allocation method.


### 4.2.2. Normal Nonlinear Allocation Method

In nonlinear allocation, other than calculating∆TE and ∆δTVN, the outputs of thrust vector engine system TEandδTVN are calculated directly through solving nonlinear equations. This method is suitable for all maneuver flights and will not cause approximation error. The total time consumption remains rational and applicable for onboard use though the calculation time is longer than the incremental linear allocation method. The flow chart of the allocation process is exhibited in Figure 6.Figure 6
Second-layer normal nonlinear allocation method.


## 4.2.1. Incremental Linear Allocation Method

The incremental linear allocation method allocates∆TE ∆δTVN and ∆δAS according to ∆x4 and ∆x2t. The incremental daisy chaining is adopted to allocate ∆x4s and ∆x4t based on (51). In this method, the priority of ∆x4s outstrips ∆x4t, and thus, the thrust vector moment is adopted only when the aerodynamic control moment gets saturated.
(51)Δx4=Δx4s+Δx4t.The increment of aerodynamic control surface is calculated by(52)ΔδAS=Δx4sq¯Sb1CAS.The reference value of aerodynamic control surface can be expressed as(53)δASref=δAS.f+ΔδAS.On that basis, introduce (49) into the incremental form, as shown in (54), in which the incremental multiplication is considered a high-order small quantity and can be ignored.
(54)Δx2tΔx4t=GEΔTEΔδTVN,(55)GE=cosδLcosδR0−TLsinδL−TRsinδR000sinδF00TFcosδFsinδLsinδR−cosδFTLcosδLTRcosδRTFsinδF−sinδLdytsinδRdyt0−TLcosδLdytTRcosδRdyt0sinδLdxtsinδRdxtcosδFdxfTLcosδLdxtTRcosδRdxt−TFsinδFdxfcosδLdyt−cosδRdytsinδFdxf−TLsinδLdytTRsinδRdytTFcosδFdxf.The thrust vector engine model is transformed from a nonlinear system into a linear system, and∆TE and ∆δTVN can be calculated directly by matrix inversion through the incremental linear allocation, as shown in (56). This method can evidently expedite allocation without iterative solution of nonlinear equations. However, like (14), when UAV is in high maneuver flight, the increment of control variables will be large, and the neglection of incremental multiplication will cause approximation errors. In this regard, this method is appropriate for minor or moderate maneuverable flights.
(56)ΔTEΔδTVN=GE−1Δx2tΔx4t.The reference thrust vector engine actuators are calculated by (57).
(57)TEref=TE.f+ΔTE,δTVNref=δTVN.f+ΔδTVN.The range ofTL, TR, TF and δL, δR, δF is limited on the basis of power system, as shown in (58). Eventually, the control signals are output to the aerodynamic control surfaces and power system.
(58)Ti=Timax,ifTi≥Timax,Ti=Timin,ifTi≤Timin,δi=minδitn−1+δ̇imaxΔt,δimax,ifδitn≥minδitn−1+δ̇imaxΔt,δimax,δi=maxδitn−1−δ̇imaxΔt,δimin,ifδitn≤maxδitn−1−δ̇imaxΔt,δimin.The flow chart of the incremental linear allocation is presented in Figure5.Figure 5
Second-layer incremental linear allocation method.


## 4.2.2. Normal Nonlinear Allocation Method

In nonlinear allocation, other than calculating∆TE and ∆δTVN, the outputs of thrust vector engine system TEandδTVN are calculated directly through solving nonlinear equations. This method is suitable for all maneuver flights and will not cause approximation error. The total time consumption remains rational and applicable for onboard use though the calculation time is longer than the incremental linear allocation method. The flow chart of the allocation process is exhibited in Figure 6.Figure 6
Second-layer normal nonlinear allocation method.


## 5. Simulation

An example is designed in this section to test the control method which is discussed in the previous section. The target trajectory is designed as follows:(59)vx=0.5t,0≤t<20,10,20≤t<60,10+t−602,60≤t<80,vy=0,0≤t<20,7sin2πt−2040,20≤t<60,0,60≤t<80,vz=−2sinπt20,0≤t<20,0,20≤t<80.The initial state of UAV is static. From 0 to 20 seconds, the vehicle is climbing up and accelerating to 10 m/s in groundx-axis. From 20 to 60 seconds, the UAV is in transition flight and maneuvering laterally, showing UAV’s flight ability in low-altitude complex flight conditions, including cities and mountainous areas. After 60 seconds, the vehicle makes uniformly accelerated flight in ground x-axis and transforms from transition flight to cruise.The example involves three different flight states, i.e., VTOL, transition flight, and cruise. These states can adequately show UAV’s longitudinal and lateral maneuverability required for an urban flight vehicle. The controller designed in this paper does not need to switch as flight states change. The trajectory can be tracked merely through adjusting the weights in the controller. As the simulation results prove, the control method can solve nonlinear, nonaffine, and coupled control problems effectively and allocate redundant control variables appropriately.To show the effectiveness of the two-layer cascaded control allocation for different flight states and mission requirements, two flight strategies are designed. In the first strategy, the UAV will track the trajectory with small attitude angle magnitude and fluctuation. In the second strategy, the vehicle is required to track the target trajectory through employing small direct force control and primarily adopting attitude control. For all strategies, the incremental linear allocation method is adopted in the second-layer control allocation, and 0.01 s is adopted as the calculation step.


### 5.1. Strategy 1: Trajectory Tracking

For Strategy1, the initial values of the weights are presented in (60). To avoid the errors caused by the small magnitudes of weighted control variables, the weights have been expanded by 100 times.
(60)W1=diag4.4122.9022.900.111.080.05,W2=diag9.5619.1319.130.240.240.24.The weight structure is generated and regulated all through the trajectory tracking, and at last, the weight changes into (61).
(61)W1=diag5.7830.2030.200.090.830.04,W2=diag6.4912.9512.950.190.190.19.Additionally, the weight change would be better consecutive for different flight states and mission requirements, since inconsecutive weight change will cause the saltation of flight states. The trajectory tracking result is presented in Figure7.Figure 7
Trajectory tracking results.The flight states of Strategy1 are exhibited in Figure 8.Figure 8
Flight states of Strategy1.
(a)
Three-axis velocity
(b)
Kinematic angle
(c)
Wind angle
(d)
Body angle
(e)
Control moments
(f)
Vectored thrust in body axis
(g)
Engine thrust
(h)
Vectoring nozzle deflection
(i)
Control surface deflection
(j)
Percentage of aerodynamic control momentsThe wind angle and body angle of the UAV are presented in Figures8(c) and 8(d). As indicated, the angle of attack α is negative in climbing phase, and α increases with the increasing of velocity and pitch angle. The UAV first keeps the minus pitch angle, and the body z-axis thrust is adopted to assist x-axis acceleration. With the rise of lift and decline of climb rate, the pitch angle (θ) and vectoring nozzle deflection angle increase, and more body x-axis thrust is adopted to accelerate the vehicle in ground x-axis. These maneuvering strategies are calculated by the controller automatically according to current flight states and weight structure.Figures8(i) and 8(j) show the changes of control surface deflections and aerodynamic control moments, respectively. In Figure 8(i), the black curve (delta-a) represents the aileron deflection angle (δa), and the blue curve (delta-e) represents the elevator deflection angle (δe). In Figure 8(j), the black curve (PTC-roll) represents the percentage of aerodynamic roll control moment in total roll control moments, and the blue curve (PTC-pitch) represents the percentage of aerodynamic pitch control moment in total pitch control moments. As the daisy chaining method is adopted in the second layer control allocation, the thrust vector control moments will be used only when aerodynamic control moments get saturated. Also, the effect of aerodynamic control surface is affected by the magnitude of wind angle. With large magnitude of angle of attack, the effect of control surface is less and vice versa. Therefore, in the takeoff phase, the control surfaces are in maximum deflection angle (±32°). From 0 to 20 seconds, with the increase of velocity, the percentage of aerodynamic control moments increases, and control surfaces are no longer saturated.From 20 to 60 seconds, the UAV is in transition flight and making lateral maneuvering, the flight velocity is 10 m/s, andα ranges from 0 degrees to 17 degrees. In Figure 8(f), Tz is negative as the body x-axis is pointing down, indicating that UAV employs vectored thrust to compensate the lift. The change tendencies of Tz and α are similar with the absolute value of velocity. As indicated, the larger lift is generated with the rise of velocity and smaller Tz is required to compensate the lift. The change of engine thrust and nozzle deflection angle after control allocation are presented in Figures 8(g) and 8(h). In Figures 8(i) and 8(j), the aerodynamic control moments get saturated again when the large control moments are required, and their usage percentage decreases as the thrust vector control moments are added in. During lateral maneuver, roll angle ϕ (Figure 8(d)) and sideslip angle β (Figure 8(c)) are changed jointly, meaning that STT and BTT controls are adopted in the UAV simultaneously. Additionally, different values of TL and TR change similarly with β. As indicated, the UAV uses engine differential thrust to realize the STT control. The change of these flight states is corresponded to kinematic azimuth angle χ, as presented in Figure 8(b). Overall, when χ is increasing, ϕ, β, and nc keep positive; when χ is decreasing, ϕ, β, and nc keep negative, as presented in Figures 8(b), 8(c), 8(d), and 8(e).Whent > 60s, with the increase of velocity, α and Tz decrease, the usage percentage of aerodynamic roll control moment rises to 100%, and the usage of aerodynamic pitch control moments increases. The UAV is transformed from transition flight to cruise.


### 5.2. Strategy 2: Trajectory Tracking

For Strategy2, the initial values of the weights are presented in
(62)W1=diag2.2822.7822.780.131.280.12,W1=diag9.9619.9219.920.280.280.28.With the generation and regulation of weight structure during trajectory tracking, at the end of simulation, the weight changes into (63).
(63)W1=diag1.8132.5332.530.080.810.08,W1=diag6.3212.6512.650.180.180.18.The flight states of Strategy2 are presented in Figure 9.Figure 9
Flight states of Strategy2.
(a)
Wind angle
(b)
Body angle
(c)
Control moments
(d)
Vectored thrust in body axis
(e)
Engine thrust
(f)
Vectoring nozzle deflectionDue to different weight selections, the tracking maneuver in Strategy2 differs from that in Strategy 1 evidently. In Strategy 2, α and θ increase quickly (Figures 9(a) and 9(b)) during takeoff, and the vehicle climbs and accelerates in forward with the thrust in body x-axis. In the lateral maneuver process, the ranges of attitude angle and control moments (Figures 9(a) and 9(c)) in Strategy 2 are larger than the one in Strategy 1 evidently. However, the vectored thrust TZ and power system actuators TF, δL, and δR adopted in Strategy 2 (Figures 9(d), 9(e), and 9(f)) are less than those in Strategy 1 (Figures 8(f), 8(g), and 8(h)). Compared with Strategy 1, Strategy 2 adopts less direct force control, and the maneuver is dominated by attitude control. In the meantime, similar with Strategy 1, BTT and STT controls are both adopted in Strategy 2. In Figure 9(c), from 60 seconds to 80 seconds, the saltation of control moment mc exists, which is caused by the aerodynamic interference between the tandem wings. During this period, α decreases from 36 degrees to 5 degrees. When α is around 18 degrees, the first saltation of pitch moment coefficient Cm is caused by the different stall angles between the front wing and the rear wing. When α is around 9 degrees, the wake flow of the front wing beats on the rear wing directly, which leads to the second saltation of Cm. As indicated, in the INDI control system, the control moment mc is used to counteract the aerodynamic disturbances and keep the stability of flight attitude.Comparing the change ofα in Figures 8(c) and 9(a), we can find that when UAV makes lateral maneuver, α in Strategy 1 increases with the increase of V, while Strategy 2 is on the contrary. This is because the average α approaching to 12 degrees is small in Strategy 1 and does not reach the stall angle. With the rise of velocity, the increment of α can largely increase lift and reduce TZ. This maneuver strategy based on weight selection is an optimal result, which makes the objective function get minor results. In Strategy 2, the average value of α is around 34 degrees, which exceeds the stall angle evidently. In this condition, increasing α will result in lift loss. Accordingly, with the increase of velocity, more lift is generated, and reducing α is the best solution for optimal control allocation.The application of the INDI control method has a requirement for the frequency of the control system. High control frequency is especially required for fiercer maneuver. In INDI control, the multiplication of increment caused by coupled control variables counts as a high-order small quantity and is neglected, which will cause large errors under fierce maneuver and large sampling time. Also, in the discrete differentiator, which is used to calculate the derivative of control vectors, the larger sampling time will result in less derivation accuracy. These errors can be restrained through reducing the sampling time, which is equivalent to increase the control frequency. As the simulation results indicate, UAV can track the target trajectory accurately. The application of the INDI method and two-layer cascaded optimal control allocation in a controller of 100 Hz can satisfy the need of maneuver in this case. Furthermore, the simulation validates the effectiveness of the control method developed in this paper.


## 5.1. Strategy 1: Trajectory Tracking

For Strategy1, the initial values of the weights are presented in (60). To avoid the errors caused by the small magnitudes of weighted control variables, the weights have been expanded by 100 times.
(60)W1=diag4.4122.9022.900.111.080.05,W2=diag9.5619.1319.130.240.240.24.The weight structure is generated and regulated all through the trajectory tracking, and at last, the weight changes into (61).
(61)W1=diag5.7830.2030.200.090.830.04,W2=diag6.4912.9512.950.190.190.19.Additionally, the weight change would be better consecutive for different flight states and mission requirements, since inconsecutive weight change will cause the saltation of flight states. The trajectory tracking result is presented in Figure7.Figure 7
Trajectory tracking results.The flight states of Strategy1 are exhibited in Figure 8.Figure 8
Flight states of Strategy1.
(a)
Three-axis velocity
(b)
Kinematic angle
(c)
Wind angle
(d)
Body angle
(e)
Control moments
(f)
Vectored thrust in body axis
(g)
Engine thrust
(h)
Vectoring nozzle deflection
(i)
Control surface deflection
(j)
Percentage of aerodynamic control momentsThe wind angle and body angle of the UAV are presented in Figures8(c) and 8(d). As indicated, the angle of attack α is negative in climbing phase, and α increases with the increasing of velocity and pitch angle. The UAV first keeps the minus pitch angle, and the body z-axis thrust is adopted to assist x-axis acceleration. With the rise of lift and decline of climb rate, the pitch angle (θ) and vectoring nozzle deflection angle increase, and more body x-axis thrust is adopted to accelerate the vehicle in ground x-axis. These maneuvering strategies are calculated by the controller automatically according to current flight states and weight structure.Figures8(i) and 8(j) show the changes of control surface deflections and aerodynamic control moments, respectively. In Figure 8(i), the black curve (delta-a) represents the aileron deflection angle (δa), and the blue curve (delta-e) represents the elevator deflection angle (δe). In Figure 8(j), the black curve (PTC-roll) represents the percentage of aerodynamic roll control moment in total roll control moments, and the blue curve (PTC-pitch) represents the percentage of aerodynamic pitch control moment in total pitch control moments. As the daisy chaining method is adopted in the second layer control allocation, the thrust vector control moments will be used only when aerodynamic control moments get saturated. Also, the effect of aerodynamic control surface is affected by the magnitude of wind angle. With large magnitude of angle of attack, the effect of control surface is less and vice versa. Therefore, in the takeoff phase, the control surfaces are in maximum deflection angle (±32°). From 0 to 20 seconds, with the increase of velocity, the percentage of aerodynamic control moments increases, and control surfaces are no longer saturated.From 20 to 60 seconds, the UAV is in transition flight and making lateral maneuvering, the flight velocity is 10 m/s, andα ranges from 0 degrees to 17 degrees. In Figure 8(f), Tz is negative as the body x-axis is pointing down, indicating that UAV employs vectored thrust to compensate the lift. The change tendencies of Tz and α are similar with the absolute value of velocity. As indicated, the larger lift is generated with the rise of velocity and smaller Tz is required to compensate the lift. The change of engine thrust and nozzle deflection angle after control allocation are presented in Figures 8(g) and 8(h). In Figures 8(i) and 8(j), the aerodynamic control moments get saturated again when the large control moments are required, and their usage percentage decreases as the thrust vector control moments are added in. During lateral maneuver, roll angle ϕ (Figure 8(d)) and sideslip angle β (Figure 8(c)) are changed jointly, meaning that STT and BTT controls are adopted in the UAV simultaneously. Additionally, different values of TL and TR change similarly with β. As indicated, the UAV uses engine differential thrust to realize the STT control. The change of these flight states is corresponded to kinematic azimuth angle χ, as presented in Figure 8(b). Overall, when χ is increasing, ϕ, β, and nc keep positive; when χ is decreasing, ϕ, β, and nc keep negative, as presented in Figures 8(b), 8(c), 8(d), and 8(e).Whent > 60s, with the increase of velocity, α and Tz decrease, the usage percentage of aerodynamic roll control moment rises to 100%, and the usage of aerodynamic pitch control moments increases. The UAV is transformed from transition flight to cruise.


## 5.2. Strategy 2: Trajectory Tracking

For Strategy2, the initial values of the weights are presented in
(62)W1=diag2.2822.7822.780.131.280.12,W1=diag9.9619.9219.920.280.280.28.With the generation and regulation of weight structure during trajectory tracking, at the end of simulation, the weight changes into (63).
(63)W1=diag1.8132.5332.530.080.810.08,W1=diag6.3212.6512.650.180.180.18.The flight states of Strategy2 are presented in Figure 9.Figure 9
Flight states of Strategy2.
(a)
Wind angle
(b)
Body angle
(c)
Control moments
(d)
Vectored thrust in body axis
(e)
Engine thrust
(f)
Vectoring nozzle deflectionDue to different weight selections, the tracking maneuver in Strategy2 differs from that in Strategy 1 evidently. In Strategy 2, α and θ increase quickly (Figures 9(a) and 9(b)) during takeoff, and the vehicle climbs and accelerates in forward with the thrust in body x-axis. In the lateral maneuver process, the ranges of attitude angle and control moments (Figures 9(a) and 9(c)) in Strategy 2 are larger than the one in Strategy 1 evidently. However, the vectored thrust TZ and power system actuators TF, δL, and δR adopted in Strategy 2 (Figures 9(d), 9(e), and 9(f)) are less than those in Strategy 1 (Figures 8(f), 8(g), and 8(h)). Compared with Strategy 1, Strategy 2 adopts less direct force control, and the maneuver is dominated by attitude control. In the meantime, similar with Strategy 1, BTT and STT controls are both adopted in Strategy 2. In Figure 9(c), from 60 seconds to 80 seconds, the saltation of control moment mc exists, which is caused by the aerodynamic interference between the tandem wings. During this period, α decreases from 36 degrees to 5 degrees. When α is around 18 degrees, the first saltation of pitch moment coefficient Cm is caused by the different stall angles between the front wing and the rear wing. When α is around 9 degrees, the wake flow of the front wing beats on the rear wing directly, which leads to the second saltation of Cm. As indicated, in the INDI control system, the control moment mc is used to counteract the aerodynamic disturbances and keep the stability of flight attitude.Comparing the change ofα in Figures 8(c) and 9(a), we can find that when UAV makes lateral maneuver, α in Strategy 1 increases with the increase of V, while Strategy 2 is on the contrary. This is because the average α approaching to 12 degrees is small in Strategy 1 and does not reach the stall angle. With the rise of velocity, the increment of α can largely increase lift and reduce TZ. This maneuver strategy based on weight selection is an optimal result, which makes the objective function get minor results. In Strategy 2, the average value of α is around 34 degrees, which exceeds the stall angle evidently. In this condition, increasing α will result in lift loss. Accordingly, with the increase of velocity, more lift is generated, and reducing α is the best solution for optimal control allocation.The application of the INDI control method has a requirement for the frequency of the control system. High control frequency is especially required for fiercer maneuver. In INDI control, the multiplication of increment caused by coupled control variables counts as a high-order small quantity and is neglected, which will cause large errors under fierce maneuver and large sampling time. Also, in the discrete differentiator, which is used to calculate the derivative of control vectors, the larger sampling time will result in less derivation accuracy. These errors can be restrained through reducing the sampling time, which is equivalent to increase the control frequency. As the simulation results indicate, UAV can track the target trajectory accurately. The application of the INDI method and two-layer cascaded optimal control allocation in a controller of 100 Hz can satisfy the need of maneuver in this case. Furthermore, the simulation validates the effectiveness of the control method developed in this paper.


## 6. Conclusions

A trajectory control system is proposed in this paper to solve the nonlinear, nonaffine, redundant, and coupled control problems existing in a UAV, which adopts direct force control technique. The research results are concluded below: (1) as the simulation indicates, the INDI controller is applicable for different flight modes, and the UAV with the INDI controller exhibits continuous maneuverability in transition flight. (2) The two-layer cascaded optimal control allocation is designed based on INDI control law; with rationally designed constraints and cascaded solving strategy, this method transforms the nonlinear coupled allocation problem into an incremental linear-equality constraint optimization problem and largely improves the allocation speed. (3) The dynamic weight strategy can adjust the weights according to different flight states and mission requirements, which ensures the rationality of control allocation. In the future, the control system presented in this paper will be incorporated with the controller of a prototype. A flight test is required to be further conducted to verify the effectiveness of the proposed control system.



---

*Source: 6315856-2018-11-27.xml*