# A Multiobjective Incremental Control Allocation Strategy for Tailless Aircraft

Authors: Maoyu Su; Jianbo Hu; Yingyang Wang; Zihou He; Jiping Cong; Linxiao Han
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

To address the control allocation problem caused by the redundant arrangement of control surfaces with nonlinear effectiveness for tailless aircraft, a novel multiobjective incremental control allocation (MICA) strategy is proposed. Firstly, the incremental nonlinear control allocation (INCA) method together with the active set quadratic programming algorithm is adopted to precisely allocate the virtual control commands. Secondly, a series of normalized objective functions in the form of increment are designed. Combining these functions by means of linear weighted sum, an incremental multiobjective function is constructed. Then, an improved nondominated sorting genetic algorithm (INSGA) is introduced to offline determine a set of weights that best meets the requirements of each flight phase. In this way, the dependence on subjective experience is minimized based on the theory of Pareto optimal. Meanwhile, the huge computational burden that the intelligent optimization algorithm brings can also be avoided. Finally, combined with the nonlinear dynamic inversion (NDI) control method, a closed-loop validation for the effectiveness of this control allocation strategy is carried out on the tailless aircraft model.


---

## Body

## 1. Introduction

With the advantages such as strong stealth and high lift-to-drag ratio [1], the tailless aircraft is a highly competitive alternative for the next-generation fighters. Due to the cancellation of vertical tail, radar cross section (RCS) is kept at a relatively low level, but the horizontal control effectiveness and lateral stability of tailless aircraft are reduced. So the tailless aircraft is usually equipped with innovative control surfaces such as all-moving wing tips and spoiler-slot deflectors besides conventional control surfaces [2]. However, this design also brings a series of challenges to flight control. First, the control surface system of tailless aircraft is highly redundant, and the deflection of each control surface will generate moments in multiple directions instead of one certain direction. So it is necessary to effectively coordinate the deflection of control surfaces to accurately and quickly generate the moments required for attitude control within the limits of deflection positions and rates, achieving the purpose of control allocation [3]. Second, the control surfaces of tailless aircraft are closely arrayed, so the aerodynamic coupling between them is serious [4]. Besides, there is often a nonlinear relationship between the single control surface deflection and the moment it generates, which further increases the difficulty of control allocation. Third, there is always more than one combination of control surface deflections to meet the requirements of flight control. The performance of different deflection combinations varies greatly, such as total energy cost of surface deflections, drag and lift. Therefore, it is necessary to optimize multiple objectives through the reasonable deflection of control surfaces according to the requirements of flight phases and tasks [5].In order to cope with the above challenges, scholars have carried out a series of relevant studies. The research on nonlinear control allocation used to deal with nonlinear control effectiveness and cross coupling effects started relatively early and developed relatively mature. It can be roughly divided into traditional methods and intelligent methods. The former is based on piecewise linearization [6], nonlinear programming [7], or other methods to overcome the allocation difficulties caused by nonlinear control effectiveness, and offset the cross-coupling effects under the guidance of compensation method [8] and robust control principle [9]. Such methods have good real-time performance, but the allocation accuracy is often limited by the operation step size, polynomial order, and other factors. The latter is based on modern intelligent algorithms such as ant colony algorithm and differential evolution algorithm [10]. Reinforcement learning [11] is also applied on nonlinear control allocation. These methods take advantage of the powerful nonlinear optimization capabilities of intelligent algorithms and neural networks to process nonlinear control effectiveness and cross-coupling effects. The high quality of solution can be ensured in most cases, but the large amount of calculation is an inevitable obstacle to accomplish its real-time application under current calculation conditions. Additionally, the relatively poor reliability of the algorithms leads to doubts about its feasibility in the field of aerospace [12].The multiobjective control allocation strategy can be roughly divided into two categories according to the idea of multiobjective optimization. The first type of strategy is based on the decomposition idea. The multibjective optimization problem is transformed into a single-objective optimization problem through methods such as linear weighted sum andε-constraint (also known as main objective method). Sun et al. [13] constructed a comprehensive objective function in the form of linear weighted sum, and used analytic hierarchy process (AHP) method to determine weights for each single objective according to the performance requirements of different flight phases. Zhang et al. [14] proposed a hybrid multiobjective control allocation strategy for the reentry phase of reusable launch vehicle, similarly using AHP method to determine the weight values in the objective function. But the weights are calculated for a series of chosen feature points rather than flight phase, and the real-time adjustment of weight values is achieved through linear interpolation. Aiming at compound high-speed rotorcraft in multimode conversion process, Zheng et al. [15] presented a hybrid multiobjective control allocation strategy. Preference matrix was introduced to determine the weights of each single objective, and an adaptive particle swarm optimization (PSO) algorithm was proposed to calculate the deflection of control surfaces with nonlinear control effectiveness; Park et al. [16] designed a multiobjective control allocation method based on evolutionary game theory, but the method is not applicable for nonlinear tailless aircraft model. Such strategies are relatively simple and easy to understand, but the determination of weights is highly subjective, and it is difficult to obtain Pareto optimal solutions under nonconvex conditions [17].The second category of the strategy is based on the Pareto optimal theory, and uses multiobjective evolutionary algorithms (MOEAs) to search for the Pareto front and Pareto optimal solutions. There are many excellent MOEAs, according to difference on the mechanism of evolution; they can be divided into 3 categories [18–21]: (1) MOEAs based on the Pareto dominance relationship; (2) MOEAs based on indicator function; and (3) MOEAs based on decomposition. The details of the 3 categories of MOEAs are concluded and compared in Table 1.Table 1
An overview of MOEAs.
Basic ideas(1) Using the Pareto fitness allocation strategy to find all the nondominated individuals from the current population(2) Using the performance evaluation indicator function to guide the search process and the solution selection process(3) Decomposing the complex multiobjective problem into several single objective subproblems or multiple simple multiobjective problemsApplication scenarios(1) Multiobjective optimization problems with 2 or 3 objectives(2) Multiobjective optimization problems with irregular Pareto front(3) Many-objective optimization problems (MaOPs)Typical algorithms(1) NSGA-II (nondominated sorting genetic algorithm-II); SPEA2 (strength Pareto evolutionary Algorithm2); PESA-II (Pareto envelop-based selection algorithm-II)(2) IBEA (indicator-based evolutionary algorithm); HypE (Hypervolume based evolutionary algorithm); SMS-EMOA (S-metric selection based evolutionary multiobjective algorithm)(3) MOEA/D (multiobjective evolutionary algorithm based on decomposition); REVA (reference vector guided evolutionary algorithm); NSGA-III (nondominated sorting genetic algorithm-III)However, in the field of multiobjective control allocation, there are few researches on using the above MOEAs to solve this problem. Taking allocation efficiency and deflection amount as optimization goals, Chen et al. [22] presented an improved NSGA to find multiple sets of Pareto optimal solutions. However, the method can only be applied to aircraft models with linear control effectiveness, and the deflection position and rate limits of control surfaces were not considered. Zhao et al. [23] designed a comprehensive multiobjective PSO algorithm and an improved hierarchical multiobjective PSO algorithm to simultaneously solve the nonlinearity of control effectiveness and multiobjective optimization problems. The advantages and disadvantages of the two methods are compared. Zheng et al. [24] proposed a hybrid multiobjective control allocation method based on indicator function based MOEAs—multiobjective PSO with adaptive probability guidance, which can comprehensively balance the optimization precision and calculation speed. To conclude, such strategies do not rely on subjective experience and can be optimized under any nonconvex conditions, but they are computationally complex. So it is necessary to take further researches on better applying MOEAs on the multiobjective control allocation problem.In this paper, a novel multiobjective nonlinear control allocation strategy is proposed. The strategy is aimed at realizing the precise allocation of control commands and the comprehensive optimization of the energy cost that control surface deflections cause, flight drag, and lift for the tailless aircraft in multiphase flight. The main contributions of this paper and the novelty of the proposed strategy are summarized as follows:(1)
A closed-loop incremental nonlinear control allocation method with the active set quadratic programming solver is presented, which is effective in precisely allocating the virtual control command against the adverse effects caused by nonlinear control effectiveness. What is more, the high precision of the proposed method is theoretically verified by proving its convergence, while other similar references [25–27] that introduce nonlinear control allocation methods usually ignore this essential step(2)
A series of objective functions are constructed to accurately characterize the performances related to control surface deflection. Then a comprehensive objective function in the form of linear weighted sum is further constructed. A critic difference of it from other comprehensive objective functions, such as the function mentioned in [13], is that each single objective function in this paper is normalized before constructing the comprehensive one, so as to facilitate the next process of weight determination(3)
A novel weight determination method based on MOEAs is proposed. Different from the commonly used expert knowledge based methods, such as the method based on AHP, the proposed method is subject to the Pareto optimal theory and takes the requirement priority for different objectives into consideration, which reduces the influence of human factor on the process of optimal weight value searching and better ensures the objectivity of the processThe remainder of this paper is organized as follows: Section2 introduces the dynamics model of the tailless aircraft and describes the multiobjective nonlinear control allocation problem; Section 3 clarifies the incremental nonlinear control allocation method [28]. Then several objective functions in the incremental form and a comprehensive multiobjective function in the form of linear weighted sum are constructed. The weight determination method of the multiobjective function based on the improved nondominated sorting genetic algorithm (INSGA) is proposed. Section 4 proves the convergence of the tailless aircraft flight control system that includes the designed multiobjective incremental control allocation (MICA) module. In Section 5, the effectiveness of the strategy is verified in the attitude control of the tailless aircraft. Section 6 presents the conclusion of this paper.


## 2. Tailless Aircraft Control Allocation Problem Description

### 2.1. Mathematical Model of Tailless Aircraft

This paper takes the Innovative Control Effectors (ICE) tailless aircraft designed by Lockheed Martin as the research object. The tailless aircraft is equipped with the multiaxis thrust vectoring and 11 aerodynamic control surfaces shown in Figure1 [29]. It is obvious that the number of actual control input far exceeds that of control moments, and some control surfaces such as spoiler slot deflector and elevon, inboard leading edge flap, and outboard leading edge flap, are closely aligned.Figure 1
ICE aircraft control surfaces configuration.Table2 introduces the deflection characteristics of each control surface.where
(1)H1s=18×100s+18s+100,H2s=40×100s+40s+100.Table 2
ICE aircraft control surfaces characteristics.
Control surfacePosition limit (°)Rate limit (°/s)DynamicsLeft/rightInboard leading edge flap(ILEF)(0,40)40H1sLeft/rightOutboard leading edge flap(OLEF)(40,40)40H1sLeft/rightAll moving wing tip(AMT)(0,60)150H2sLeft/right Elevon(ELE)(-30,30)150H2sLeft/rightSpoiler slot deflector(SSD)(0,60)150H2sPitch flap(PF)(-30,30)150H2sPitch/yawMultiaxis thrust vectoring(MTV)(-15,15)150H2sThe rotational motion model and the center-of-mass dynamics model of the ICE tailless aircraft in the airflow coordinate system are the following [30];
(2)α̇=−pcosαtanβ+q−rsinαtanβ+cosμcosβgvcosγ−FXsinαmvcosβ+FZcosαmvcosβ,β̇=psinα−rcosα+gvsinμcosγ−FXcosαsinβmv+FYcosβmv−FZsinαsinβmv,μ̇=pcosα+rsinαcosβ−gvcosγcosμtanβ+FYtanγcosβcosμmv+tanγsinμ+tanβmvFXsinα−FZcosα−tanγcosμsinβmvFXcosα+FZsinα,(3)ṗ=IzL+IzxN+IzxIx−Iy+Izpq+IyIz−Izx2−Iz2qrIxIz−Izx2,q̇=M+Izxr2−p2+Iz−IxprIy,ṙ=IzxL+IxN+Ix2−IxIy+Izx2pq+IzxIy−Iz−IxqrIxIz−Izx2,where p, q, and r are the rolling, pitching, and yawing angular velocities; L, M, N, and Ix, Iy, Iz are the three-axis control moment and the corresponding inertia moment of the aircraft in the body coordinate system, respectively; Izx is the product of inertia; FX, FY and FZ are the components of the resultant force that includes engine thrust and aerodynamic force along the body axes; α, β, μ are, respectively, the aerodynamic angle of attack, the aerodynamic angle of sideslip, and the rolling angle around the velocity vector axis; γ is the track inclination in the track coordinate system.For the tailless aircraft, there is a complex nonlinear mapping relationship between the deflection of the control surface and its generated moment. Take the pitching moment coefficient described by the wind tunnel test data of the ICE aircraft as an example:(4)Cm=Mq¯Sc¯=Cm_baseα,Ma+q¯c¯2VCm_q¯α,Ma+Cm_betaα,β,Ma+fflxercmalt,MaCm_leleα,δlssd,δlele,Ma+Cm_releα,δrssd,δrele,Ma+Cm_lilefα,β,δlilef+Cm_rilefα,β,δrilef+Cm_lolefα,β,δlilef,δlolef,Ma+Cm_rolefα,β,δrilef,δrolef,Ma+Cm_lamt−lolefα,δlolef,δlamt+Cm_ramt−rolefα,δrolef,δramt+Cm_lamt−leleα,δlele,δlamt+Cm_ramt−releα,δrele,δramt+Cm_pfα,β,δrssd,δlssd,δpf,Ma+Cm_lamtα,β,δlamt+Cm_ramtα,β,δramt+fssd−pfMaCm_lssdα,β,δlssd+Cm_rssdα,β,δrssd,where Ma is the Mach number, q¯ is the dynamic pressure, S is the wing reference area, and c¯ is the wing mean geometric chord. For the meanings of the other symbols, please refer to [29]. From Equation (4) and Figure 1, it can be seen that the control effectiveness varies due to the close arrayment between the control surfaces such as the SSD and ELE. The control effectiveness will be affected by the deflection of adjacent control surfaces to varying degrees, and the cross-coupling effects are too obvious to be neglected. Figure 2 directly reflects the huge difference and nonlinear characteristics of the control effectiveness of each control surface on the ICE aircraft. For example, the control effectiveness curve of SSD is of monotonic nonlinearity and the control effectiveness curves of Cl−δamt and Cm−δamt are of nonmonotonic nonlinearity.Figure 2
Control effectiveness curve of control surfaces.
(a)
Rolling moment generation effectiveness(b)
Pitching moment generation effectiveness(c)
Yawing moment generation effectiveness


### 2.2. Nonlinear Dynamic Inversion Control Law Design

According to the separation principle of time scale [31], that is, the controlled variables can be divided into fast variables and slow variables based on the their response speed and different control laws can be independently designed for the two kinds of variables, the flight control system is divided into inner and outer loops, which have been shown in Equation (2) and Equation (3), marked as
(5)ẋ1=f1x1+g1x1x2,ẋ2=f2x2+g2x2τ,y=x1,x2T,where τ=LMNT is the virtual control command, x1=αβμT is regarded as the slow variables and x2=pqrT is regarded as the fast variables.The outer loop and inner loop dynamic inversion control laws are designed as(6)x2=g1−1x1ω1x1c−x1−f1x1,τ=g2−1x2ω2x2c−x2−f2x2,where ω1 and ω2 are diagonal matrices, and their diagonal elements are all positive real numbers.


### 2.3. Description of Multiobjective Nonlinear Control Allocation Problem

In order to coordinate the deflection of control surfaces more effectively and reduce the design difficulty of the flight control system, the flight control system of the tailless aircraft with redundant control surfaces usually adopts the modular design scheme shown in Figure3.Figure 3
Modular flight control system structure of tailless aircraft.First, the flight control law calculates the virtual control commandτc∈Rm according to the state feedback, in this paper, τc=CldCmdCndT, m=3, Cld=Ld/q¯Sb, Cmd=Md/q¯Sc¯, Cnd=Nd/q¯Sb are, respectively, the desired roll, pitch, and yaw moment coefficient, and b is the wing span. Then, based on the nonlinear control effectiveness model
(7)τ=hx,δ,and the position and rate limits of the control surface deflection
(8)δmin≤δ≤δmax,δ̇≤δ̇max,where δmax and δmin are the upper and lower limits of the deflection position of the control surface, respectively, δ̇max is the maximum deflection rate of the control surface, the control allocation module runs a certain control allocation algorithm to calculate the control surface deflection command δc. The actual deflection of control surfaces is δ, generating the virtual control input τ=ClCmCnT and finally changing the flight state x. In many studies, in order to simplify subsequent processing, Equation (7) is linearly approximated as
(9)τ≈Bxδ≈Bδ,where Bx∈Rm×p is called the control effectiveness matrix, p is the number of control surface. Then the control allocation problem can be solved through traditional methods such as direct allocation, redistribution pseudo inverse, and linear programming. However, as mentioned in Section 2.1, the control surfaces on the ICE tailless aircraft have complex nonlinear control characteristics, so Equation (7) cannot be approximated to the form of Equation (9), and it is impossible to use the above methods to calculate the control surface deflections. The characteristics greatly increase the difficulty of mapping from virtual control commands to control surface deflection commands.In addition, control allocation can also be regarded as a constrained optimization problem, and the optimization objective is(10)minJ=t−tc2=hx,d−tc2.In order to ensure the uniqueness of the solution, one or several secondary objectives related to the deflection of control surfaces are usually added. It is constructed as a multiobjective optimization function in the form of(11)minJ=FJ1δ,J2δ,⋯,Jnδ,where n is the number of the secondary objective. It should be noted that in this multiobjective optimization problem, the minimum allocation error must always be the primary optimization objective to ensure the stability of the control system.However, due to the complex aerodynamic characteristics of tailless aircraft, some secondary objectives are difficult to precisely depict, such as drag and lift. Due to the diversity of tasks performed by the tailless aircraft, it is necessary to dynamically adjust the priority of each secondary objective in the multiobjective optimization function. Besides, the interaction relationship between secondary objectives further increases the difficulty of this multiobjective optimization problem.


## 2.1. Mathematical Model of Tailless Aircraft

This paper takes the Innovative Control Effectors (ICE) tailless aircraft designed by Lockheed Martin as the research object. The tailless aircraft is equipped with the multiaxis thrust vectoring and 11 aerodynamic control surfaces shown in Figure1 [29]. It is obvious that the number of actual control input far exceeds that of control moments, and some control surfaces such as spoiler slot deflector and elevon, inboard leading edge flap, and outboard leading edge flap, are closely aligned.Figure 1
ICE aircraft control surfaces configuration.Table2 introduces the deflection characteristics of each control surface.where
(1)H1s=18×100s+18s+100,H2s=40×100s+40s+100.Table 2
ICE aircraft control surfaces characteristics.
Control surfacePosition limit (°)Rate limit (°/s)DynamicsLeft/rightInboard leading edge flap(ILEF)(0,40)40H1sLeft/rightOutboard leading edge flap(OLEF)(40,40)40H1sLeft/rightAll moving wing tip(AMT)(0,60)150H2sLeft/right Elevon(ELE)(-30,30)150H2sLeft/rightSpoiler slot deflector(SSD)(0,60)150H2sPitch flap(PF)(-30,30)150H2sPitch/yawMultiaxis thrust vectoring(MTV)(-15,15)150H2sThe rotational motion model and the center-of-mass dynamics model of the ICE tailless aircraft in the airflow coordinate system are the following [30];
(2)α̇=−pcosαtanβ+q−rsinαtanβ+cosμcosβgvcosγ−FXsinαmvcosβ+FZcosαmvcosβ,β̇=psinα−rcosα+gvsinμcosγ−FXcosαsinβmv+FYcosβmv−FZsinαsinβmv,μ̇=pcosα+rsinαcosβ−gvcosγcosμtanβ+FYtanγcosβcosμmv+tanγsinμ+tanβmvFXsinα−FZcosα−tanγcosμsinβmvFXcosα+FZsinα,(3)ṗ=IzL+IzxN+IzxIx−Iy+Izpq+IyIz−Izx2−Iz2qrIxIz−Izx2,q̇=M+Izxr2−p2+Iz−IxprIy,ṙ=IzxL+IxN+Ix2−IxIy+Izx2pq+IzxIy−Iz−IxqrIxIz−Izx2,where p, q, and r are the rolling, pitching, and yawing angular velocities; L, M, N, and Ix, Iy, Iz are the three-axis control moment and the corresponding inertia moment of the aircraft in the body coordinate system, respectively; Izx is the product of inertia; FX, FY and FZ are the components of the resultant force that includes engine thrust and aerodynamic force along the body axes; α, β, μ are, respectively, the aerodynamic angle of attack, the aerodynamic angle of sideslip, and the rolling angle around the velocity vector axis; γ is the track inclination in the track coordinate system.For the tailless aircraft, there is a complex nonlinear mapping relationship between the deflection of the control surface and its generated moment. Take the pitching moment coefficient described by the wind tunnel test data of the ICE aircraft as an example:(4)Cm=Mq¯Sc¯=Cm_baseα,Ma+q¯c¯2VCm_q¯α,Ma+Cm_betaα,β,Ma+fflxercmalt,MaCm_leleα,δlssd,δlele,Ma+Cm_releα,δrssd,δrele,Ma+Cm_lilefα,β,δlilef+Cm_rilefα,β,δrilef+Cm_lolefα,β,δlilef,δlolef,Ma+Cm_rolefα,β,δrilef,δrolef,Ma+Cm_lamt−lolefα,δlolef,δlamt+Cm_ramt−rolefα,δrolef,δramt+Cm_lamt−leleα,δlele,δlamt+Cm_ramt−releα,δrele,δramt+Cm_pfα,β,δrssd,δlssd,δpf,Ma+Cm_lamtα,β,δlamt+Cm_ramtα,β,δramt+fssd−pfMaCm_lssdα,β,δlssd+Cm_rssdα,β,δrssd,where Ma is the Mach number, q¯ is the dynamic pressure, S is the wing reference area, and c¯ is the wing mean geometric chord. For the meanings of the other symbols, please refer to [29]. From Equation (4) and Figure 1, it can be seen that the control effectiveness varies due to the close arrayment between the control surfaces such as the SSD and ELE. The control effectiveness will be affected by the deflection of adjacent control surfaces to varying degrees, and the cross-coupling effects are too obvious to be neglected. Figure 2 directly reflects the huge difference and nonlinear characteristics of the control effectiveness of each control surface on the ICE aircraft. For example, the control effectiveness curve of SSD is of monotonic nonlinearity and the control effectiveness curves of Cl−δamt and Cm−δamt are of nonmonotonic nonlinearity.Figure 2
Control effectiveness curve of control surfaces.
(a)
Rolling moment generation effectiveness(b)
Pitching moment generation effectiveness(c)
Yawing moment generation effectiveness


## 2.2. Nonlinear Dynamic Inversion Control Law Design

According to the separation principle of time scale [31], that is, the controlled variables can be divided into fast variables and slow variables based on the their response speed and different control laws can be independently designed for the two kinds of variables, the flight control system is divided into inner and outer loops, which have been shown in Equation (2) and Equation (3), marked as
(5)ẋ1=f1x1+g1x1x2,ẋ2=f2x2+g2x2τ,y=x1,x2T,where τ=LMNT is the virtual control command, x1=αβμT is regarded as the slow variables and x2=pqrT is regarded as the fast variables.The outer loop and inner loop dynamic inversion control laws are designed as(6)x2=g1−1x1ω1x1c−x1−f1x1,τ=g2−1x2ω2x2c−x2−f2x2,where ω1 and ω2 are diagonal matrices, and their diagonal elements are all positive real numbers.


## 2.3. Description of Multiobjective Nonlinear Control Allocation Problem

In order to coordinate the deflection of control surfaces more effectively and reduce the design difficulty of the flight control system, the flight control system of the tailless aircraft with redundant control surfaces usually adopts the modular design scheme shown in Figure3.Figure 3
Modular flight control system structure of tailless aircraft.First, the flight control law calculates the virtual control commandτc∈Rm according to the state feedback, in this paper, τc=CldCmdCndT, m=3, Cld=Ld/q¯Sb, Cmd=Md/q¯Sc¯, Cnd=Nd/q¯Sb are, respectively, the desired roll, pitch, and yaw moment coefficient, and b is the wing span. Then, based on the nonlinear control effectiveness model
(7)τ=hx,δ,and the position and rate limits of the control surface deflection
(8)δmin≤δ≤δmax,δ̇≤δ̇max,where δmax and δmin are the upper and lower limits of the deflection position of the control surface, respectively, δ̇max is the maximum deflection rate of the control surface, the control allocation module runs a certain control allocation algorithm to calculate the control surface deflection command δc. The actual deflection of control surfaces is δ, generating the virtual control input τ=ClCmCnT and finally changing the flight state x. In many studies, in order to simplify subsequent processing, Equation (7) is linearly approximated as
(9)τ≈Bxδ≈Bδ,where Bx∈Rm×p is called the control effectiveness matrix, p is the number of control surface. Then the control allocation problem can be solved through traditional methods such as direct allocation, redistribution pseudo inverse, and linear programming. However, as mentioned in Section 2.1, the control surfaces on the ICE tailless aircraft have complex nonlinear control characteristics, so Equation (7) cannot be approximated to the form of Equation (9), and it is impossible to use the above methods to calculate the control surface deflections. The characteristics greatly increase the difficulty of mapping from virtual control commands to control surface deflection commands.In addition, control allocation can also be regarded as a constrained optimization problem, and the optimization objective is(10)minJ=t−tc2=hx,d−tc2.In order to ensure the uniqueness of the solution, one or several secondary objectives related to the deflection of control surfaces are usually added. It is constructed as a multiobjective optimization function in the form of(11)minJ=FJ1δ,J2δ,⋯,Jnδ,where n is the number of the secondary objective. It should be noted that in this multiobjective optimization problem, the minimum allocation error must always be the primary optimization objective to ensure the stability of the control system.However, due to the complex aerodynamic characteristics of tailless aircraft, some secondary objectives are difficult to precisely depict, such as drag and lift. Due to the diversity of tasks performed by the tailless aircraft, it is necessary to dynamically adjust the priority of each secondary objective in the multiobjective optimization function. Besides, the interaction relationship between secondary objectives further increases the difficulty of this multiobjective optimization problem.


## 3. Multiobjective Incremental Control Allocation Strategy

### 3.1. Incremental Nonlinear Control Allocation

In order to solve the problem that linear control allocation methods cannot solve due to the nonlinearity of control effectiveness and the cross-coupling of control surfaces, this paper adopts the incremental nonlinear control allocation (INCA) method to deal with it.The first-order Taylor expansion of Equation (7) can be carried out at x0,δ0.
(12)τ=τ0+Δτ=hx0,δ0+∂h∂xδ=δ0x=x0x−x0+∂h∂δδ=δ0x=x0δ−δ0,where x0, δ0, and τ0 are the flight state, the deflection position of control surfaces, and the virtual control input at last step, respectively.Since the sampling frequency of the flight control computer is usually around 100 Hz, so the sampling interval is small enough, according to the separation principle of time scale, the flight statex is a slow variable relative to the deflection δ of control surfaces, and the state x can be regarded as a constant value in the duration of about 0.01 s. Therefore, the expression of Δτ in Equation (12) can be simplified into the form of multiplying the Jacobian matrix by the difference between the adjacent two-step control surface deflections.
(13)Δτ=∂∂δhx,δδ=δ0x=x0Δδ=∇δhx0,δ0Δδ,where
(14)∇δhx,δ=∂h1x,δ∂δ1∂h1x,δ∂δ2···∂h1x,δ∂δp∂h2x,δ∂δ1∂h2x,δ∂δ2···∂h2x,δ∂δp············∂hmx,δ∂δ1∂hmx,δ∂δ2···∂hmx,δ∂δp,and ∇δhx0,δ0 is also called the local control effectiveness matrix at the deflection position δ0, indicating the matrix changes as the control surfaces deflect rather than keep constant like the matrix B in Equation (9). Each partial derivative in Equation (14) can be approximated by the difference method,
(15)∂h∂δi≈hx0,δ1,δ2,⋯,δi+Δ,⋯−hx0,δ1,δ2,⋯,δi−Δ,⋯2Δ,where Δ is a very small amount, which is taken as 0.01° in this paper. Meanwhile, the position and rate limits of the control surface deflection can also be rewritten in the incremental form:
(16)Δδmin≤Δδ≤Δδmax,Δδmax=minδ̇maxΔts,δmax−δ0,Δδmin=max−δ̇maxΔts,δmin−δ0,where Δts is the sampling step. The minimum allocation error objective function of the control allocation problem in the incremental form is
(17)minJ0=∇δhx,δΔδ+τ0−τc2,Δδmin≤Δδ≤Δδmax.On this basis, the redistribution pseudo inverse and mathematical programming algorithms can be used to calculate the control surfaces deflection increment. Then the total deflection can be obtained by adding it to the deflection of last step. Figure4 summarizes the procedure of INCA method:Figure 4
Flow chart of incremental nonlinear control allocation method.According to the description, it can be seen that INCA has the following advantages:(1)
High precision. Similar to the mechanism of the feedback control system that eliminates the error between output and input, the closed-loop structure enables the allocation error to be eliminated through allocating the difference signal between virtual control command and virtual control input generated in last operation cycle. In Section4, the convergence of this method will be proved, so as to verify its high precision advantage theoretically(2)
Good real-time performance. INCA is based on the discrete operation characteristics of the flight control computer. It obtains the local control effectiveness matrix by calculating the Jacobian matrix in real time. Then, the nonlinear control allocation problem is transformed into a linear one in each operation cycle, making the efficient linear control allocation methods still applicable(3)
Control surface limits can be effectively handled. The deflection position and rate limits are unified into the incremental form shown in Equation (16), reducing the difficulty of handling the control surface limits, especially the deflection rate limits


### 3.2. Construction of Secondary Objectives in Incremental Form

After completing the incremental form reconstruction of the nonlinear control allocation problem, the multiobjective optimization problem needs to be further considered. This paper plans to comprehensively optimize four objectives, including allocation error, total surface deflections, flight drag, and lift produced by control surface deflection. First, the incremental form objective function is constructed for each secondary objective, respectively.


#### 3.2.1. Minimum Surface Deflection

The minimum surface deflection is the most commonly used secondary objective in the control allocation problem, and together with the minimum allocation error, it constitutes a hybrid optimization objective. Achieving the minimum surface deflection objective is of practical significance in three aspects: (1) it can reduce the control energy consumption caused by the deflection of the control surface; (2) it can prevent the control surface from reaching the saturation position prematurely and improve the control allocation efficiency [32]; (3) the deflection of control surface effects the RCS of aircraft [33], and relatively little deflection keeps the RCS as low as possible. The expression of this objective in the incremental form is
(18)minJ1=Δδ+δ02,Δδmin≤Δδ≤Δδmax.


#### 3.2.2. Minimum Drag

Due to the complex aerodynamic characteristics of tailless aircraft, its drag effectiveness is difficult to describe. For example, under the flight conditions ofMa=1.2 and h=15000ft, the drag effectiveness curve of each control surface is shown in Figure 5:Figure 5
Control surface drag effectiveness curve of tailless aircraft.Some studies choose the minimum surface deflection objective to approximate the minimum drag objective. However, due to the difference in the characteristics of different control surfaces, there is still large room to reduce the drag generated by the deflection of control surfaces.It is known that the dragD, lift L, and lateral force Y of the aircraft can be obtained from the following equations [34]
(19)DLY=−cosαcosβ−cosαsinβ−sinαsinαcosβsinαsinβ−cosα−sinβcosβ0FX−TFYFZ.In this paper, referring to [12], the drag coefficient is written in the incremental form.
(20)CDx,δ=CD0x,δ0+∇δCDx,δ0Δδ,(21)∇δCDx,δ0=∂CD∂δ1x,δ0⋯⋯∂CD∂δpx,δ0.Similar to Equation (15), each partial derivative in Equation (21) can also be approximated by the difference form.
(22)∂CD∂δi=CDx,δ1,δ2,⋯,δi+Δ,⋯−CDx,δ1,δ2,⋯,δi−Δ,⋯2Δ.The minimum drag objective is then constructed as follows:(23)minJ2=∇δCDx,δ0Δδ+CD02,Δδmin≤Δδ≤Δδmax.


#### 3.2.3. Maximum Lift

Similar to the drag effectiveness, the lift effectiveness of the tailless vehicle are also complex, which is depicted in Figure6.Figure 6
Lift effectiveness curve of control surface of tailless aircraft.The maximum lift objective is constructed as(24)minJ3=∇δCLx,δ0Δδ+CL0−λ2,Δδmin≤Δδ≤Δδmax,where λ is a positive real number that satisfies
(25)∇δCLx,δ0Δδ+CL0−λ<0,∀Δδ∈Δδmin,Δδmax.


#### 3.2.4. Comprehensive Objective

On this basis, aiming to comprehensively optimize the above objectives, and according to the performance requirements of the tailless aircraft in different flight phases, a comprehensive objective function is further constructed. To facilitate the solution, the comprehensive objective function is constructed in the form of linear weighted sum as follows:(26)J=kJ0+∑i=13εiJiHowever, due to the great difference in the order of magnitude of each optimization objective, it is difficult to determine the weight in the comprehensive objective function in the form of linear weighted sum. So it is necessary to normalize each single objective in the comprehensive function in the following form as follows:(27)J0∗=J0/τ0−τc2max,J1∗=J1/δ2max,J2∗=J2/CDmax,J3∗=J3/λ,J∗=kJ0∗+∑i=13εiJi∗,where ·2max represents the maximum that the 2-norm of a variable can reach within its upper and lower limits. In order to use the active set quadratic programming algorithm to solve the problem, the normalized comprehensive objective function can be written in the form of weighted least squares:
(28)J∗=FΔδ+G2=kτ0−τc2max∇δhx,δ0ε1δ2maxε2CDmax∇δCDx,δ0ε3λ∇δCLx,δ0Δδ+kτ0−τc2maxτ0−τcε1δ2maxδε2CDmaxCD0ε3λCL0−λ2.When the weightsk, ε1, ε2, and ε3 are given, the increment of surface deflection can be directly obtained by running the active set quadratic programming algorithm, and then the total amount of surface deflection can be calculated. It can be seen that the determination of weights is of great importance and needs to be achieved according to the flight phase and task requirements. Nevertheless, each objective is not independent of each other, that is, the optimization of one objective might obstruct or promote the optimization of other objectives, as qualitatively shown in Figure 7.Figure 7
Qualitative relationship between some secondary objectives.To be specific, when the weightk of the minimum error objective is fixed, the increase of the weight of any other second objective will reduce the consideration of the minimum error when calculating the deflection of the control surface through active set quadratic programming. In addition, because the drag and lift are partially derived from the deflection of the control surface, for the objectives of minimum drag and maximum lift established above, when the weight of the objective of minimum total deflection of the control surface increases, the drag and lift will remain at a relatively low level. Therefore, the minimum control surface deflection and minimum drag objective promote each other, while the objective of minimum control surface deflection and the objective of maximum lift restrain each other. To conclude, there is a complex unproportionate relationship between the weight values and flight performance, which is shown in Figure 7.The weight determination method based on AHP never takes the interaction between secondary objectives into consideration, and highly relies on expert experience, which is of great subjectivity. So the optimization result is often not the Pareto optimal solution that best meets the requirements. Therefore, according to the Pareto dominance relationship and the priority relationship of each optimization objective in different flight phases, this paper adopts an intelligent weight determination method of more objectivity based on the improved NSGA algorithm to achieve multiobjective control allocation.


### 3.3. Multiobjective Weight Determination Based on Improved NSGA

Owing to the powerful nonlinear optimization capabilities, metaheuristic optimization algorithms are suitable for the determination of weights to achieve the optimization ofJ1J2J3. However, due to the unbearable computational burden, it is difficult for the current flight control computer to support real-time solution based on meta-heuristic optimization algorithms. For the purpose of fostering strengths and avoiding weaknesses, this paper adopts MOEAs to offline determine the weight of each secondary objective for different flight phases and tasks. Furthermore, the problem of weight determination in the multiobjective function has the following characteristics: (1) this is a multiobjective optimization problem with 3 objectives; (2) the true Pareto front is unknown; and (3) the preference for each objective changes with the flight phase. Therefore, according to the information of basic ideas and application scenarios in Table 1, NSGA-II [35] is selected among so many MOEAs mentioned in Section 1, and, in consideration of different preference for each objective, it is improved by modifying the expression of crowding distance to optimize the weight in Equation (28).Above all, some basic concepts need to be clarified:Definition 1 (Pareto dominance).
In a multiobjective minimum optimization problem, assume thatε1 and ε2 are two sets of feasible solution within the limits. It is called that ε1 Pareto dominates ε2 if and only if ∀i∈1,2,⋯,n, Jiε1≤Jiε2, and ∃i∈1,2,⋯,n, Jiε1<Jiε2, which is noted as ε1≺ε2.Definition 2 (Pareto optimal solution).
In the setΩ formed by all feasible solutions within the limits, if a feasible solution ε∗ is not dominated by any other feasible solution, then ε∗ is called the Pareto optimal solution of the multiobjective optimization problem.Definition 3 (Pareto optimal solution set).
There is often more than one Pareto optimal solution, and the set of all Pareto optimal solutions is called the Pareto optimal solution set, denoted asP.Definition 4 (Pareto front).
The set of function values of all Pareto optimal solutions in the objective space is called the Pareto front, denoted asPf=J1ε∗,J2ε∗,⋯,Jnε∗ε∗∈P.The procedure to determine the weight of each objective in the multiobjective control allocation problem through improved NSGA-II (INSGA) is as follows:Step 5.
Set the population sizeN and the max iteration imax, generate the initial parent population P0 through random initialization. Then perform nondominated sorting to generate several classification subsets F1, F2,···, Fl, where F1 is the set of all nondominated individuals in P0, F2 is the set of all nondominated individuals in P0−F1, and so on. Also define a partial order relationship between two individuals.
For two individualsε1∈Fi and ε2∈Fj in the population P0, the relationship between their partial orders Poε1 and Poε2 is
(29)Poε1>Poε2ifi<jorDε1>Dε2andi=j.
D· represents the crowding distance of the individual, and Poε1>Poε2 represents that the partial order of the individual ε1 is better than the partial order of the individual ε2. ε1 will be preferentially selected in the selection operation.
The standard NSGA-II algorithm defines the crowding distance as(30)Dεi=∑j=1nJjεi+1−Jjεi−1,where εi+1 and εi−1 are the two individuals most adjacent to εi, and the crowding distance of the boundary point is set to ∞. However, because the priority w of each objective is different under different circumstances, though sorting by the crowding distance shown in Equation (30) can maintain the diversity of the population, it is not conducive to ensure that the high-quality individuals better meeting the requirements of different phases are retained. So the crowding distance is improved to the following form:
(31)Dεi=∑j=1nJjεi+1−Jjεi−1wjJjεi.
In the flight of the tailless aircraft, the prioritieswj of different objectives at each phase are shown in Table 3, where 1 represents the relatively highest priority, and 3 represents the relatively lowest priority.Step 6.
According to the nondominated sorting results and the partial order relationship, the parent populationPii=0,1,2,⋯,imax−1 is selected in the method of binary championship, crossover, and mutation are operated to create the son population Si in the same size N. The process of crossover operation is
(32)ai1j=rand⋅ai2j+1−rand⋅ai1j,ai2j=rand⋅ai1j+1−rand⋅ai2j,where ai1j and ai2j are the jth genes of individuals numbered i1 and i2 in the parent population Pi, i1,i2=1,2,⋯,p and i1≠i2, j=1,2,⋯,q, and rand is a random number in the interval (0,1). The process of mutation operation is
(33)aij=aij+aij−amax⋅rand⋅1−iterimax2,rand>0.5,aij−amin−aij⋅rand⋅1−iterimax2,rand≤0.5,where iter is the current number of iterations, and the influence of the mutation operation will be gradually weakened with the increase of the number of iterations. The design plays a role in improving the robustness and global optimality of the algorithm.Step 7.
Merge the parent population and the son population into a populationQi in the size of 2 N, and perform nondominated sorting on it.Step 8.
Determine the partial order relationship between individuals in the critical layer classification subset according to the crowding distance, where the definition of critical layer classification subset is: in several classification subsets ofQi, the number nFk of individuals contained in subset Fk satisfies
(34)∑i=1k−1nFi<N,∑i=1knFi≥N,then Fk is called the critical classification subset of Qi.Step 9.
According to the nondominated sorting result and the partial order relationship, selectN individuals to form the parent population Pi+1 in next generation. The specific process is shown in Figure 8.Step 10.
To obtain the final populationPimax and its corresponding Pareto front, repeat step 6 to step 9 until the max iteration is reached.Step 11.
In the application scenario of this paper, it is necessary to select an optimal individualPimaxk from Pimax as the weights in Equation (28). According to the priority relationship of each optimization objective corresponding to each flight phase predetermined in Table 2, first, for Ji that has the highest priority, select the top a% of the optimal individuals from Pimax to form a loose set A. Then, for the next highest priority objective Jj, select the top b% of the optimal individuals from the loose set A to form a new loose set B. Finally, for Jk with the lowest priority objective, select the optimal individual in B. Take the individual as the weight of each subobjective in the objective function shown in Equation (28).Table 3
Priority of each secondary objective.
Secondary objectiveClimbingCruiseManeuverLandingMinimum deflection3212Minimum drag2133Maximum lift1321Figure 8
New population generation mechanism in INSGA.The overall flow of the algorithm is summarized in Figure9:Figure 9
Flow chart of INSGA.


### 3.4. Tailless Aircraft Control System with MICA Module

According to the above description of nonlinear dynamic inversion, INCA method and multiobjective weight determination method based on INSGA, a complete strategy of objective incremental control allocation (MICA) for tailless aircraft is formed. Figure10 summarizes the tailless aircraft flight control system adopting the MICA strategy.Figure 10
Block diagram of tailless aircraft control system with MICA module.


## 3.1. Incremental Nonlinear Control Allocation

In order to solve the problem that linear control allocation methods cannot solve due to the nonlinearity of control effectiveness and the cross-coupling of control surfaces, this paper adopts the incremental nonlinear control allocation (INCA) method to deal with it.The first-order Taylor expansion of Equation (7) can be carried out at x0,δ0.
(12)τ=τ0+Δτ=hx0,δ0+∂h∂xδ=δ0x=x0x−x0+∂h∂δδ=δ0x=x0δ−δ0,where x0, δ0, and τ0 are the flight state, the deflection position of control surfaces, and the virtual control input at last step, respectively.Since the sampling frequency of the flight control computer is usually around 100 Hz, so the sampling interval is small enough, according to the separation principle of time scale, the flight statex is a slow variable relative to the deflection δ of control surfaces, and the state x can be regarded as a constant value in the duration of about 0.01 s. Therefore, the expression of Δτ in Equation (12) can be simplified into the form of multiplying the Jacobian matrix by the difference between the adjacent two-step control surface deflections.
(13)Δτ=∂∂δhx,δδ=δ0x=x0Δδ=∇δhx0,δ0Δδ,where
(14)∇δhx,δ=∂h1x,δ∂δ1∂h1x,δ∂δ2···∂h1x,δ∂δp∂h2x,δ∂δ1∂h2x,δ∂δ2···∂h2x,δ∂δp············∂hmx,δ∂δ1∂hmx,δ∂δ2···∂hmx,δ∂δp,and ∇δhx0,δ0 is also called the local control effectiveness matrix at the deflection position δ0, indicating the matrix changes as the control surfaces deflect rather than keep constant like the matrix B in Equation (9). Each partial derivative in Equation (14) can be approximated by the difference method,
(15)∂h∂δi≈hx0,δ1,δ2,⋯,δi+Δ,⋯−hx0,δ1,δ2,⋯,δi−Δ,⋯2Δ,where Δ is a very small amount, which is taken as 0.01° in this paper. Meanwhile, the position and rate limits of the control surface deflection can also be rewritten in the incremental form:
(16)Δδmin≤Δδ≤Δδmax,Δδmax=minδ̇maxΔts,δmax−δ0,Δδmin=max−δ̇maxΔts,δmin−δ0,where Δts is the sampling step. The minimum allocation error objective function of the control allocation problem in the incremental form is
(17)minJ0=∇δhx,δΔδ+τ0−τc2,Δδmin≤Δδ≤Δδmax.On this basis, the redistribution pseudo inverse and mathematical programming algorithms can be used to calculate the control surfaces deflection increment. Then the total deflection can be obtained by adding it to the deflection of last step. Figure4 summarizes the procedure of INCA method:Figure 4
Flow chart of incremental nonlinear control allocation method.According to the description, it can be seen that INCA has the following advantages:(1)
High precision. Similar to the mechanism of the feedback control system that eliminates the error between output and input, the closed-loop structure enables the allocation error to be eliminated through allocating the difference signal between virtual control command and virtual control input generated in last operation cycle. In Section4, the convergence of this method will be proved, so as to verify its high precision advantage theoretically(2)
Good real-time performance. INCA is based on the discrete operation characteristics of the flight control computer. It obtains the local control effectiveness matrix by calculating the Jacobian matrix in real time. Then, the nonlinear control allocation problem is transformed into a linear one in each operation cycle, making the efficient linear control allocation methods still applicable(3)
Control surface limits can be effectively handled. The deflection position and rate limits are unified into the incremental form shown in Equation (16), reducing the difficulty of handling the control surface limits, especially the deflection rate limits


## 3.2. Construction of Secondary Objectives in Incremental Form

After completing the incremental form reconstruction of the nonlinear control allocation problem, the multiobjective optimization problem needs to be further considered. This paper plans to comprehensively optimize four objectives, including allocation error, total surface deflections, flight drag, and lift produced by control surface deflection. First, the incremental form objective function is constructed for each secondary objective, respectively.


### 3.2.1. Minimum Surface Deflection

The minimum surface deflection is the most commonly used secondary objective in the control allocation problem, and together with the minimum allocation error, it constitutes a hybrid optimization objective. Achieving the minimum surface deflection objective is of practical significance in three aspects: (1) it can reduce the control energy consumption caused by the deflection of the control surface; (2) it can prevent the control surface from reaching the saturation position prematurely and improve the control allocation efficiency [32]; (3) the deflection of control surface effects the RCS of aircraft [33], and relatively little deflection keeps the RCS as low as possible. The expression of this objective in the incremental form is
(18)minJ1=Δδ+δ02,Δδmin≤Δδ≤Δδmax.


### 3.2.2. Minimum Drag

Due to the complex aerodynamic characteristics of tailless aircraft, its drag effectiveness is difficult to describe. For example, under the flight conditions ofMa=1.2 and h=15000ft, the drag effectiveness curve of each control surface is shown in Figure 5:Figure 5
Control surface drag effectiveness curve of tailless aircraft.Some studies choose the minimum surface deflection objective to approximate the minimum drag objective. However, due to the difference in the characteristics of different control surfaces, there is still large room to reduce the drag generated by the deflection of control surfaces.It is known that the dragD, lift L, and lateral force Y of the aircraft can be obtained from the following equations [34]
(19)DLY=−cosαcosβ−cosαsinβ−sinαsinαcosβsinαsinβ−cosα−sinβcosβ0FX−TFYFZ.In this paper, referring to [12], the drag coefficient is written in the incremental form.
(20)CDx,δ=CD0x,δ0+∇δCDx,δ0Δδ,(21)∇δCDx,δ0=∂CD∂δ1x,δ0⋯⋯∂CD∂δpx,δ0.Similar to Equation (15), each partial derivative in Equation (21) can also be approximated by the difference form.
(22)∂CD∂δi=CDx,δ1,δ2,⋯,δi+Δ,⋯−CDx,δ1,δ2,⋯,δi−Δ,⋯2Δ.The minimum drag objective is then constructed as follows:(23)minJ2=∇δCDx,δ0Δδ+CD02,Δδmin≤Δδ≤Δδmax.


### 3.2.3. Maximum Lift

Similar to the drag effectiveness, the lift effectiveness of the tailless vehicle are also complex, which is depicted in Figure6.Figure 6
Lift effectiveness curve of control surface of tailless aircraft.The maximum lift objective is constructed as(24)minJ3=∇δCLx,δ0Δδ+CL0−λ2,Δδmin≤Δδ≤Δδmax,where λ is a positive real number that satisfies
(25)∇δCLx,δ0Δδ+CL0−λ<0,∀Δδ∈Δδmin,Δδmax.


### 3.2.4. Comprehensive Objective

On this basis, aiming to comprehensively optimize the above objectives, and according to the performance requirements of the tailless aircraft in different flight phases, a comprehensive objective function is further constructed. To facilitate the solution, the comprehensive objective function is constructed in the form of linear weighted sum as follows:(26)J=kJ0+∑i=13εiJiHowever, due to the great difference in the order of magnitude of each optimization objective, it is difficult to determine the weight in the comprehensive objective function in the form of linear weighted sum. So it is necessary to normalize each single objective in the comprehensive function in the following form as follows:(27)J0∗=J0/τ0−τc2max,J1∗=J1/δ2max,J2∗=J2/CDmax,J3∗=J3/λ,J∗=kJ0∗+∑i=13εiJi∗,where ·2max represents the maximum that the 2-norm of a variable can reach within its upper and lower limits. In order to use the active set quadratic programming algorithm to solve the problem, the normalized comprehensive objective function can be written in the form of weighted least squares:
(28)J∗=FΔδ+G2=kτ0−τc2max∇δhx,δ0ε1δ2maxε2CDmax∇δCDx,δ0ε3λ∇δCLx,δ0Δδ+kτ0−τc2maxτ0−τcε1δ2maxδε2CDmaxCD0ε3λCL0−λ2.When the weightsk, ε1, ε2, and ε3 are given, the increment of surface deflection can be directly obtained by running the active set quadratic programming algorithm, and then the total amount of surface deflection can be calculated. It can be seen that the determination of weights is of great importance and needs to be achieved according to the flight phase and task requirements. Nevertheless, each objective is not independent of each other, that is, the optimization of one objective might obstruct or promote the optimization of other objectives, as qualitatively shown in Figure 7.Figure 7
Qualitative relationship between some secondary objectives.To be specific, when the weightk of the minimum error objective is fixed, the increase of the weight of any other second objective will reduce the consideration of the minimum error when calculating the deflection of the control surface through active set quadratic programming. In addition, because the drag and lift are partially derived from the deflection of the control surface, for the objectives of minimum drag and maximum lift established above, when the weight of the objective of minimum total deflection of the control surface increases, the drag and lift will remain at a relatively low level. Therefore, the minimum control surface deflection and minimum drag objective promote each other, while the objective of minimum control surface deflection and the objective of maximum lift restrain each other. To conclude, there is a complex unproportionate relationship between the weight values and flight performance, which is shown in Figure 7.The weight determination method based on AHP never takes the interaction between secondary objectives into consideration, and highly relies on expert experience, which is of great subjectivity. So the optimization result is often not the Pareto optimal solution that best meets the requirements. Therefore, according to the Pareto dominance relationship and the priority relationship of each optimization objective in different flight phases, this paper adopts an intelligent weight determination method of more objectivity based on the improved NSGA algorithm to achieve multiobjective control allocation.


## 3.2.1. Minimum Surface Deflection

The minimum surface deflection is the most commonly used secondary objective in the control allocation problem, and together with the minimum allocation error, it constitutes a hybrid optimization objective. Achieving the minimum surface deflection objective is of practical significance in three aspects: (1) it can reduce the control energy consumption caused by the deflection of the control surface; (2) it can prevent the control surface from reaching the saturation position prematurely and improve the control allocation efficiency [32]; (3) the deflection of control surface effects the RCS of aircraft [33], and relatively little deflection keeps the RCS as low as possible. The expression of this objective in the incremental form is
(18)minJ1=Δδ+δ02,Δδmin≤Δδ≤Δδmax.


## 3.2.2. Minimum Drag

Due to the complex aerodynamic characteristics of tailless aircraft, its drag effectiveness is difficult to describe. For example, under the flight conditions ofMa=1.2 and h=15000ft, the drag effectiveness curve of each control surface is shown in Figure 5:Figure 5
Control surface drag effectiveness curve of tailless aircraft.Some studies choose the minimum surface deflection objective to approximate the minimum drag objective. However, due to the difference in the characteristics of different control surfaces, there is still large room to reduce the drag generated by the deflection of control surfaces.It is known that the dragD, lift L, and lateral force Y of the aircraft can be obtained from the following equations [34]
(19)DLY=−cosαcosβ−cosαsinβ−sinαsinαcosβsinαsinβ−cosα−sinβcosβ0FX−TFYFZ.In this paper, referring to [12], the drag coefficient is written in the incremental form.
(20)CDx,δ=CD0x,δ0+∇δCDx,δ0Δδ,(21)∇δCDx,δ0=∂CD∂δ1x,δ0⋯⋯∂CD∂δpx,δ0.Similar to Equation (15), each partial derivative in Equation (21) can also be approximated by the difference form.
(22)∂CD∂δi=CDx,δ1,δ2,⋯,δi+Δ,⋯−CDx,δ1,δ2,⋯,δi−Δ,⋯2Δ.The minimum drag objective is then constructed as follows:(23)minJ2=∇δCDx,δ0Δδ+CD02,Δδmin≤Δδ≤Δδmax.


## 3.2.3. Maximum Lift

Similar to the drag effectiveness, the lift effectiveness of the tailless vehicle are also complex, which is depicted in Figure6.Figure 6
Lift effectiveness curve of control surface of tailless aircraft.The maximum lift objective is constructed as(24)minJ3=∇δCLx,δ0Δδ+CL0−λ2,Δδmin≤Δδ≤Δδmax,where λ is a positive real number that satisfies
(25)∇δCLx,δ0Δδ+CL0−λ<0,∀Δδ∈Δδmin,Δδmax.


## 3.2.4. Comprehensive Objective

On this basis, aiming to comprehensively optimize the above objectives, and according to the performance requirements of the tailless aircraft in different flight phases, a comprehensive objective function is further constructed. To facilitate the solution, the comprehensive objective function is constructed in the form of linear weighted sum as follows:(26)J=kJ0+∑i=13εiJiHowever, due to the great difference in the order of magnitude of each optimization objective, it is difficult to determine the weight in the comprehensive objective function in the form of linear weighted sum. So it is necessary to normalize each single objective in the comprehensive function in the following form as follows:(27)J0∗=J0/τ0−τc2max,J1∗=J1/δ2max,J2∗=J2/CDmax,J3∗=J3/λ,J∗=kJ0∗+∑i=13εiJi∗,where ·2max represents the maximum that the 2-norm of a variable can reach within its upper and lower limits. In order to use the active set quadratic programming algorithm to solve the problem, the normalized comprehensive objective function can be written in the form of weighted least squares:
(28)J∗=FΔδ+G2=kτ0−τc2max∇δhx,δ0ε1δ2maxε2CDmax∇δCDx,δ0ε3λ∇δCLx,δ0Δδ+kτ0−τc2maxτ0−τcε1δ2maxδε2CDmaxCD0ε3λCL0−λ2.When the weightsk, ε1, ε2, and ε3 are given, the increment of surface deflection can be directly obtained by running the active set quadratic programming algorithm, and then the total amount of surface deflection can be calculated. It can be seen that the determination of weights is of great importance and needs to be achieved according to the flight phase and task requirements. Nevertheless, each objective is not independent of each other, that is, the optimization of one objective might obstruct or promote the optimization of other objectives, as qualitatively shown in Figure 7.Figure 7
Qualitative relationship between some secondary objectives.To be specific, when the weightk of the minimum error objective is fixed, the increase of the weight of any other second objective will reduce the consideration of the minimum error when calculating the deflection of the control surface through active set quadratic programming. In addition, because the drag and lift are partially derived from the deflection of the control surface, for the objectives of minimum drag and maximum lift established above, when the weight of the objective of minimum total deflection of the control surface increases, the drag and lift will remain at a relatively low level. Therefore, the minimum control surface deflection and minimum drag objective promote each other, while the objective of minimum control surface deflection and the objective of maximum lift restrain each other. To conclude, there is a complex unproportionate relationship between the weight values and flight performance, which is shown in Figure 7.The weight determination method based on AHP never takes the interaction between secondary objectives into consideration, and highly relies on expert experience, which is of great subjectivity. So the optimization result is often not the Pareto optimal solution that best meets the requirements. Therefore, according to the Pareto dominance relationship and the priority relationship of each optimization objective in different flight phases, this paper adopts an intelligent weight determination method of more objectivity based on the improved NSGA algorithm to achieve multiobjective control allocation.


## 3.3. Multiobjective Weight Determination Based on Improved NSGA

Owing to the powerful nonlinear optimization capabilities, metaheuristic optimization algorithms are suitable for the determination of weights to achieve the optimization ofJ1J2J3. However, due to the unbearable computational burden, it is difficult for the current flight control computer to support real-time solution based on meta-heuristic optimization algorithms. For the purpose of fostering strengths and avoiding weaknesses, this paper adopts MOEAs to offline determine the weight of each secondary objective for different flight phases and tasks. Furthermore, the problem of weight determination in the multiobjective function has the following characteristics: (1) this is a multiobjective optimization problem with 3 objectives; (2) the true Pareto front is unknown; and (3) the preference for each objective changes with the flight phase. Therefore, according to the information of basic ideas and application scenarios in Table 1, NSGA-II [35] is selected among so many MOEAs mentioned in Section 1, and, in consideration of different preference for each objective, it is improved by modifying the expression of crowding distance to optimize the weight in Equation (28).Above all, some basic concepts need to be clarified:Definition 1 (Pareto dominance).
In a multiobjective minimum optimization problem, assume thatε1 and ε2 are two sets of feasible solution within the limits. It is called that ε1 Pareto dominates ε2 if and only if ∀i∈1,2,⋯,n, Jiε1≤Jiε2, and ∃i∈1,2,⋯,n, Jiε1<Jiε2, which is noted as ε1≺ε2.Definition 2 (Pareto optimal solution).
In the setΩ formed by all feasible solutions within the limits, if a feasible solution ε∗ is not dominated by any other feasible solution, then ε∗ is called the Pareto optimal solution of the multiobjective optimization problem.Definition 3 (Pareto optimal solution set).
There is often more than one Pareto optimal solution, and the set of all Pareto optimal solutions is called the Pareto optimal solution set, denoted asP.Definition 4 (Pareto front).
The set of function values of all Pareto optimal solutions in the objective space is called the Pareto front, denoted asPf=J1ε∗,J2ε∗,⋯,Jnε∗ε∗∈P.The procedure to determine the weight of each objective in the multiobjective control allocation problem through improved NSGA-II (INSGA) is as follows:Step 5.
Set the population sizeN and the max iteration imax, generate the initial parent population P0 through random initialization. Then perform nondominated sorting to generate several classification subsets F1, F2,···, Fl, where F1 is the set of all nondominated individuals in P0, F2 is the set of all nondominated individuals in P0−F1, and so on. Also define a partial order relationship between two individuals.
For two individualsε1∈Fi and ε2∈Fj in the population P0, the relationship between their partial orders Poε1 and Poε2 is
(29)Poε1>Poε2ifi<jorDε1>Dε2andi=j.
D· represents the crowding distance of the individual, and Poε1>Poε2 represents that the partial order of the individual ε1 is better than the partial order of the individual ε2. ε1 will be preferentially selected in the selection operation.
The standard NSGA-II algorithm defines the crowding distance as(30)Dεi=∑j=1nJjεi+1−Jjεi−1,where εi+1 and εi−1 are the two individuals most adjacent to εi, and the crowding distance of the boundary point is set to ∞. However, because the priority w of each objective is different under different circumstances, though sorting by the crowding distance shown in Equation (30) can maintain the diversity of the population, it is not conducive to ensure that the high-quality individuals better meeting the requirements of different phases are retained. So the crowding distance is improved to the following form:
(31)Dεi=∑j=1nJjεi+1−Jjεi−1wjJjεi.
In the flight of the tailless aircraft, the prioritieswj of different objectives at each phase are shown in Table 3, where 1 represents the relatively highest priority, and 3 represents the relatively lowest priority.Step 6.
According to the nondominated sorting results and the partial order relationship, the parent populationPii=0,1,2,⋯,imax−1 is selected in the method of binary championship, crossover, and mutation are operated to create the son population Si in the same size N. The process of crossover operation is
(32)ai1j=rand⋅ai2j+1−rand⋅ai1j,ai2j=rand⋅ai1j+1−rand⋅ai2j,where ai1j and ai2j are the jth genes of individuals numbered i1 and i2 in the parent population Pi, i1,i2=1,2,⋯,p and i1≠i2, j=1,2,⋯,q, and rand is a random number in the interval (0,1). The process of mutation operation is
(33)aij=aij+aij−amax⋅rand⋅1−iterimax2,rand>0.5,aij−amin−aij⋅rand⋅1−iterimax2,rand≤0.5,where iter is the current number of iterations, and the influence of the mutation operation will be gradually weakened with the increase of the number of iterations. The design plays a role in improving the robustness and global optimality of the algorithm.Step 7.
Merge the parent population and the son population into a populationQi in the size of 2 N, and perform nondominated sorting on it.Step 8.
Determine the partial order relationship between individuals in the critical layer classification subset according to the crowding distance, where the definition of critical layer classification subset is: in several classification subsets ofQi, the number nFk of individuals contained in subset Fk satisfies
(34)∑i=1k−1nFi<N,∑i=1knFi≥N,then Fk is called the critical classification subset of Qi.Step 9.
According to the nondominated sorting result and the partial order relationship, selectN individuals to form the parent population Pi+1 in next generation. The specific process is shown in Figure 8.Step 10.
To obtain the final populationPimax and its corresponding Pareto front, repeat step 6 to step 9 until the max iteration is reached.Step 11.
In the application scenario of this paper, it is necessary to select an optimal individualPimaxk from Pimax as the weights in Equation (28). According to the priority relationship of each optimization objective corresponding to each flight phase predetermined in Table 2, first, for Ji that has the highest priority, select the top a% of the optimal individuals from Pimax to form a loose set A. Then, for the next highest priority objective Jj, select the top b% of the optimal individuals from the loose set A to form a new loose set B. Finally, for Jk with the lowest priority objective, select the optimal individual in B. Take the individual as the weight of each subobjective in the objective function shown in Equation (28).Table 3
Priority of each secondary objective.
Secondary objectiveClimbingCruiseManeuverLandingMinimum deflection3212Minimum drag2133Maximum lift1321Figure 8
New population generation mechanism in INSGA.The overall flow of the algorithm is summarized in Figure9:Figure 9
Flow chart of INSGA.


## 3.4. Tailless Aircraft Control System with MICA Module

According to the above description of nonlinear dynamic inversion, INCA method and multiobjective weight determination method based on INSGA, a complete strategy of objective incremental control allocation (MICA) for tailless aircraft is formed. Figure10 summarizes the tailless aircraft flight control system adopting the MICA strategy.Figure 10
Block diagram of tailless aircraft control system with MICA module.


## 4. Proof of Convergence

### 4.1. Closed-Loop Convergence Proof for Incremental Nonlinear Control Allocation

According to the description of INCA in Section3.1, denoting ∇δhx0,δ0 as Bdk, the active set quadratic programming algorithm can be regarded as a pseudoinverse operation, which is expressed as
(35)Δδk=Bd+kΔτk.Then the actually generated virtual control input is(36)τk=Bakδk−1+Bd+kΔτck,where Bak is the actual control effectiveness matrix at the current moment. Assuming that between two adjacent sampling points, the change of Bak is negligible, that is, Bak=Bak−1. Denote BakBd+k as Kk, then Equation (36) can be transformed into
(37)τk=Bak−1δk−1+BakBd+kΔτck=τk−1+Kkτck−τk−1.Performz transform on Equation (37).
(38)τz=z−1τz+Kkτcz−z−1τz.Further simplify Equation (38) into
(39)τz=z−1I+Kk−1zKkτcz.The necessary and sufficient condition to ensure the stability of the closed-loop system shown in Equation (39) is that every solution of characteristic equation
(40)z−1I+Kk=0,is in the unit circle [36].According to Schur’s lemma, there are an invertible matrixP and an upper triangular matrix Q that satisfy
(41)Kk=PQP−1=Pλ1a12···a1mλ2···a2m······λmP−1.Equation (40) can be transformed into
(42)z−1PP−1+PQP−1=0⇔Pz−1I+QP−1=0⇔z−1I+Q=0.Further solve(43)zi−1+λi=0,i=1,2,⋯,m,to get
(44)zi=1−λi.Then, the necessary and sufficient condition to ensure the convergence of the closed-loop system is(45)zi=1−λi<1Setλi=xi+jyi, j is an imaginary unit, then we can get
(46)1−xi−jyi=xi−12+yi2<1⇔xi−12+yi2<1.It can be obtained that the necessary and sufficient condition for the stability of the closed-loop system is that the characteristic roots ofKk are all in the circle with the center of (1,0) and the radius of 1 in the complex plane. In other words, the closed-loop system can be converged by rationally designing Bd+k so that each characteristic root of Kk satisfies the condition of Equation (45).Applying the final value theorem ofz transform to Equation (39) can further calculate the steady-state error of the closed-loop system: let the virtual control command τck be the unit step signal, then the final value of the output virtual control variable τk is
(47)limk⟶∞τk=limz⟶1z−1τz,limk⟶∞τk=limz⟶1z−1z−1I+Kk−1zKkzz−1zz−1zz−1T,limk⟶∞τk=limz⟶1z−1I+Kk−1Kkz2z2z2T,limk⟶∞τk=111T.It can be seen that the INCA method using the active set quadratic programming solver can perform error-free allocation of virtual control commands in the form of step while ensuring the convergence of the system.


### 4.2. Convergence Proof for Nonlinear Dynamic Inversion

For the outer loop, after substituting the dynamic inversion control law into the state equation, we get(48)ẋ1=ω1x1c−x1.Since the diagonal elements ofω1 are all positive numbers, each eigenvalue of ‐ω1 contains a negative real part, and the control command x1c is bounded, so the outer loop control system is stable.For the inner loop, there will always be a certain error in the control allocation module, so the actual control input isτ=τc−e, and for the allocation error e, τ=g2−1x2ω2x2c−x2−f2x2−e is substituted into the state equation of the inner loop:
(49)ẋ2=−ω2x2+ω2x2c−g2x2e.From the closed-loop stability proof of INCA in Section4.1, it can be seen that the allocation error e is bounded and globally convergent, so we have
(50)ẋ2+ω2x2=ω2x2c−g2x2e<R+,where R+ is a positive real number. Equation (50) makes the following conditions hold:
(51)ẋ2i<0,∀x2i>R+ω2i,ẋ2i>0,∀x2i<−R+ω2i,and get
(52)x2i≤R+ω2i,This shows thatx2 is bounded and x2⟶x2c when e⟶0.


## 4.1. Closed-Loop Convergence Proof for Incremental Nonlinear Control Allocation

According to the description of INCA in Section3.1, denoting ∇δhx0,δ0 as Bdk, the active set quadratic programming algorithm can be regarded as a pseudoinverse operation, which is expressed as
(35)Δδk=Bd+kΔτk.Then the actually generated virtual control input is(36)τk=Bakδk−1+Bd+kΔτck,where Bak is the actual control effectiveness matrix at the current moment. Assuming that between two adjacent sampling points, the change of Bak is negligible, that is, Bak=Bak−1. Denote BakBd+k as Kk, then Equation (36) can be transformed into
(37)τk=Bak−1δk−1+BakBd+kΔτck=τk−1+Kkτck−τk−1.Performz transform on Equation (37).
(38)τz=z−1τz+Kkτcz−z−1τz.Further simplify Equation (38) into
(39)τz=z−1I+Kk−1zKkτcz.The necessary and sufficient condition to ensure the stability of the closed-loop system shown in Equation (39) is that every solution of characteristic equation
(40)z−1I+Kk=0,is in the unit circle [36].According to Schur’s lemma, there are an invertible matrixP and an upper triangular matrix Q that satisfy
(41)Kk=PQP−1=Pλ1a12···a1mλ2···a2m······λmP−1.Equation (40) can be transformed into
(42)z−1PP−1+PQP−1=0⇔Pz−1I+QP−1=0⇔z−1I+Q=0.Further solve(43)zi−1+λi=0,i=1,2,⋯,m,to get
(44)zi=1−λi.Then, the necessary and sufficient condition to ensure the convergence of the closed-loop system is(45)zi=1−λi<1Setλi=xi+jyi, j is an imaginary unit, then we can get
(46)1−xi−jyi=xi−12+yi2<1⇔xi−12+yi2<1.It can be obtained that the necessary and sufficient condition for the stability of the closed-loop system is that the characteristic roots ofKk are all in the circle with the center of (1,0) and the radius of 1 in the complex plane. In other words, the closed-loop system can be converged by rationally designing Bd+k so that each characteristic root of Kk satisfies the condition of Equation (45).Applying the final value theorem ofz transform to Equation (39) can further calculate the steady-state error of the closed-loop system: let the virtual control command τck be the unit step signal, then the final value of the output virtual control variable τk is
(47)limk⟶∞τk=limz⟶1z−1τz,limk⟶∞τk=limz⟶1z−1z−1I+Kk−1zKkzz−1zz−1zz−1T,limk⟶∞τk=limz⟶1z−1I+Kk−1Kkz2z2z2T,limk⟶∞τk=111T.It can be seen that the INCA method using the active set quadratic programming solver can perform error-free allocation of virtual control commands in the form of step while ensuring the convergence of the system.


## 4.2. Convergence Proof for Nonlinear Dynamic Inversion

For the outer loop, after substituting the dynamic inversion control law into the state equation, we get(48)ẋ1=ω1x1c−x1.Since the diagonal elements ofω1 are all positive numbers, each eigenvalue of ‐ω1 contains a negative real part, and the control command x1c is bounded, so the outer loop control system is stable.For the inner loop, there will always be a certain error in the control allocation module, so the actual control input isτ=τc−e, and for the allocation error e, τ=g2−1x2ω2x2c−x2−f2x2−e is substituted into the state equation of the inner loop:
(49)ẋ2=−ω2x2+ω2x2c−g2x2e.From the closed-loop stability proof of INCA in Section4.1, it can be seen that the allocation error e is bounded and globally convergent, so we have
(50)ẋ2+ω2x2=ω2x2c−g2x2e<R+,where R+ is a positive real number. Equation (50) makes the following conditions hold:
(51)ẋ2i<0,∀x2i>R+ω2i,ẋ2i>0,∀x2i<−R+ω2i,and get
(52)x2i≤R+ω2i,This shows thatx2 is bounded and x2⟶x2c when e⟶0.


## 5. Simulation and Analysis

### 5.1. Simulation Settings

The cruise phase of ICE tailless aircraft is taken as an example to conduct simulation experiments. The simulation step sizeΔts is set to 0.01 s, Ma=0.6, h=15000 ft, and the attitude tracking command is set to
(53)αdest=10lt−5lt−2.5,βdest=0,μdest=5lt−lt−1.8−lt−3.2,where lt is the unit step function.Due to the mechanism of elitist preservation in INSGA, there is no need to worry about the destruction of high-quality individuals in the crossover and mutation operations. Therefore, the crossover probabilityPc and mutation probability Pm can be set relatively high to improve the global search ability and speed up the search progress. The parameter settings of INSGA are shown in Table 4.Table 4
Algorithm parameter settings.
Population sizeMax iterationCrossover probabilityPcMutation probabilityPmRetained percentagea%Retained percentageb%504010.320%20%


### 5.2. Simulation Results and Analysis

First, run INSGA offline to solve the Pareto front of multiobjective weight and its corresponding Pareto optimal solution in the set flight state. The results are shown in Figure11.Figure 11
Pareto front and its side views.
(a)
Pareto front(b)
CL‐CD side view(c)
δ‐CD side view(d)
δ‐CL side viewFigure11 can verify the qualitative relationship between the objectives shown in Figure 7. The maximum lift and the minimum drag target are mutually restrained; within a certain range, the minimum drag and the minimum surface deflection objective are positively correlated, and the maximum lift and the minimum surface deflection are negatively correlated. According to the demand priority of each objective in the cruise phase, an optimal point (marked in red in Figure 11) is selected on the Pareto front, and its corresponding independent variable, that is, the weight of each objective, is got. Its corresponding weight ε=0.97410.00710.0188.On the condition that the model used is also the ICE aircraft and the solution algorithm used is also active set quadratic programming, the multiobjective weight of the cruise phase determined by the AHP method in Reference [13] is ε=0.54980.36810.0821. Figure 12 is a comparative analysis of the effects of the weight determined by the two methods.Figure 12
Tracking of angle commands.
(a)
Tracking of angle of attack(b)
Tracking of sideslip angle(c)
Tracking of roll angleSince the MICA method used in this paper is mainly aimed at minimizing the allocation error, the two groups of experiments with different weights can both achieve good tracking results for the attitude commands. The moment command allocation results shown in Figure13 further verify the effective processing capability of the INCA method for nonlinear control effectiveness, and the allocation errors are all within the allowable range of flight control.Figure 13
Results of moment command allocation and allocation error.
(a)
Allocation result ofCl using INCA-AHP(b)
Allocation result ofCl using INCA-INSGA(c)
Allocation result ofCm using INCA-AHP(d)
Allocation result ofCm using INCA-INSGA(e)
Allocation result ofCn using INCA-AHP(f)
Allocation result ofCn using INCA-INSGAFigure14 shows the variation of the deflection position of each control surface on the ICE aircraft during the control command tracking process. Figure 15 makes statistics on the average deflection of each control surface. It can be clearly seen that the control surface deflections differ greatly when different multiobjective weights are used. The deflection of each control surface obtained by MICA based on INSGA is generally less than 50% of that of the other method, indicating that the control energy consumption is effectively reduced. Due to the difference in control surface deflections, combined with the drag and lift generation effectiveness of each control surface of the ICE aircraft shown in Figures 5 and 6, it can be seen that due to more usage of SSDs with high drag generation effectiveness and negative lift effectiveness, the allocation results obtained by AHP-based MICA generate relatively large drag and relatively small lift. The drag and lift curves shown in Figures 16 and 17 validate this. Table 5 makes a summarized comparison on the performance of the two MICAs: the average total surface deflection of INSGA-MICA is only 31.2% of that of AHP-MICA, the average drag is 5.88% lower than that of AHP-MICA, and the average lift is 9.19% higher than that of AHP-MICA. In other words, the weights calculated by INSGA-MICA dominate the weights calculated by AHP-MICA on the performance of multiple secondary objective optimizations.Figure 14
Control surface deflection curve.Figure 15
Statistics of the average deflection of the control surface.Figure 16
Variation curve of drag coefficient.Figure 17
Variation curve of lift coefficient.Table 5
Performance comparison of AHP-MICA and INSGA-MICA.
Methoderror_Cxmaxerror_CxavgδavgCDavgCLavgAHP-MICA6.32×10−34.46×10−449.35401.206×10−20.20141.88×10−27.16×10−43.46×10−31.51×10−4INSGA-MICA6.36×10−34.22×10−415.39941.135×10−20.21992.09×10−26.90×10−43.65×10−31.74×10−4


## 5.1. Simulation Settings

The cruise phase of ICE tailless aircraft is taken as an example to conduct simulation experiments. The simulation step sizeΔts is set to 0.01 s, Ma=0.6, h=15000 ft, and the attitude tracking command is set to
(53)αdest=10lt−5lt−2.5,βdest=0,μdest=5lt−lt−1.8−lt−3.2,where lt is the unit step function.Due to the mechanism of elitist preservation in INSGA, there is no need to worry about the destruction of high-quality individuals in the crossover and mutation operations. Therefore, the crossover probabilityPc and mutation probability Pm can be set relatively high to improve the global search ability and speed up the search progress. The parameter settings of INSGA are shown in Table 4.Table 4
Algorithm parameter settings.
Population sizeMax iterationCrossover probabilityPcMutation probabilityPmRetained percentagea%Retained percentageb%504010.320%20%


## 5.2. Simulation Results and Analysis

First, run INSGA offline to solve the Pareto front of multiobjective weight and its corresponding Pareto optimal solution in the set flight state. The results are shown in Figure11.Figure 11
Pareto front and its side views.
(a)
Pareto front(b)
CL‐CD side view(c)
δ‐CD side view(d)
δ‐CL side viewFigure11 can verify the qualitative relationship between the objectives shown in Figure 7. The maximum lift and the minimum drag target are mutually restrained; within a certain range, the minimum drag and the minimum surface deflection objective are positively correlated, and the maximum lift and the minimum surface deflection are negatively correlated. According to the demand priority of each objective in the cruise phase, an optimal point (marked in red in Figure 11) is selected on the Pareto front, and its corresponding independent variable, that is, the weight of each objective, is got. Its corresponding weight ε=0.97410.00710.0188.On the condition that the model used is also the ICE aircraft and the solution algorithm used is also active set quadratic programming, the multiobjective weight of the cruise phase determined by the AHP method in Reference [13] is ε=0.54980.36810.0821. Figure 12 is a comparative analysis of the effects of the weight determined by the two methods.Figure 12
Tracking of angle commands.
(a)
Tracking of angle of attack(b)
Tracking of sideslip angle(c)
Tracking of roll angleSince the MICA method used in this paper is mainly aimed at minimizing the allocation error, the two groups of experiments with different weights can both achieve good tracking results for the attitude commands. The moment command allocation results shown in Figure13 further verify the effective processing capability of the INCA method for nonlinear control effectiveness, and the allocation errors are all within the allowable range of flight control.Figure 13
Results of moment command allocation and allocation error.
(a)
Allocation result ofCl using INCA-AHP(b)
Allocation result ofCl using INCA-INSGA(c)
Allocation result ofCm using INCA-AHP(d)
Allocation result ofCm using INCA-INSGA(e)
Allocation result ofCn using INCA-AHP(f)
Allocation result ofCn using INCA-INSGAFigure14 shows the variation of the deflection position of each control surface on the ICE aircraft during the control command tracking process. Figure 15 makes statistics on the average deflection of each control surface. It can be clearly seen that the control surface deflections differ greatly when different multiobjective weights are used. The deflection of each control surface obtained by MICA based on INSGA is generally less than 50% of that of the other method, indicating that the control energy consumption is effectively reduced. Due to the difference in control surface deflections, combined with the drag and lift generation effectiveness of each control surface of the ICE aircraft shown in Figures 5 and 6, it can be seen that due to more usage of SSDs with high drag generation effectiveness and negative lift effectiveness, the allocation results obtained by AHP-based MICA generate relatively large drag and relatively small lift. The drag and lift curves shown in Figures 16 and 17 validate this. Table 5 makes a summarized comparison on the performance of the two MICAs: the average total surface deflection of INSGA-MICA is only 31.2% of that of AHP-MICA, the average drag is 5.88% lower than that of AHP-MICA, and the average lift is 9.19% higher than that of AHP-MICA. In other words, the weights calculated by INSGA-MICA dominate the weights calculated by AHP-MICA on the performance of multiple secondary objective optimizations.Figure 14
Control surface deflection curve.Figure 15
Statistics of the average deflection of the control surface.Figure 16
Variation curve of drag coefficient.Figure 17
Variation curve of lift coefficient.Table 5
Performance comparison of AHP-MICA and INSGA-MICA.
Methoderror_Cxmaxerror_CxavgδavgCDavgCLavgAHP-MICA6.32×10−34.46×10−449.35401.206×10−20.20141.88×10−27.16×10−43.46×10−31.51×10−4INSGA-MICA6.36×10−34.22×10−415.39941.135×10−20.21992.09×10−26.90×10−43.65×10−31.74×10−4


## 6. Conclusion

A novel multiobjective control allocation strategy for tailless aircraft is designed based on the INCA method and the INSGA algorithm. The strategy mainly solves the following problems:(1)
Construct INCA framework based on piecewise linearization and error feedback, effectively overcoming the control allocation difficulties caused by the nonlinear control effectiveness and deflection limits. Furthermore, active set quadratic programming algorithm is adopted to solve the amount of deflection. The closed-loop convergence of INCA is also proved(2)
Construct an incremental linear weighted multiobjective function concerning allocation error, control surface deflections, drag, and lift. Run the INSGA to determine a set of optimal weights is selected and apply it to the constructed objective function. On the basis, Multiobjective Incremental Control Allocation strategy is completely formed(3)
Compare the advantages and disadvantages of the multiobjective weights determined by the INSGA and AHP. The simulation takes ICE tailless aircraft in the cruise phase as an example. Its results show that the proposed control allocation strategy can help generate lower the flight drag and total surface deflection, together with a higher liftThe strategy proposed in this paper comprehensively considers the accuracy and real-time requirements of the control allocation algorithm, always takes the minimum allocation error as the primary objective, and uses the multiobjective intelligent optimization algorithm to search the optimal weight offline and use it online. While exploiting the control potential of tailless aircraft, it avoids the large amount of computation and long computation time. In the next step of research, the minimum RCS objective can be further accurately characterized and included in the consideration of multiobjective optimization to meet the stealth requirements in the penetration mission of tailless aircraft. With the capability of flight computer improved, the intelligent weight determination can be operated online, so that the multiobjective weights can be adjusted more flexibly in the complex and changeable flight process and better adapt to the task environment.



---

*Source: 6515234-2022-12-16.xml*