Fuzzy-genetic Control of Quadrotor Unmanned Aerial Vehicles

This article presents a novel fuzzy identification method for dynamic modelling of quadrotor unmanned aerial vehicles. The method is based on a special parameterization of the antecedent part of fuzzy systems that results in fuzzy-partitions for antecedents. This antecedent parameter representation method of fuzzy rules ensures upholding of predefined linguistic value ordering and ensures that fuzzy-partitions remain intact throughout an unconstrained hybrid evolutionary and gradient descent based optimization process. In the equations of motion the first order derivative component is calculated based on Christoffel symbols, the derivatives of fuzzy systems are used for modelling the Coriolis effects, gyroscopic and centrifugal terms. The non-linear parameters are subjected to an initial global evolutionary optimization scheme and fine tuning with gradient descent based local search. Simulation results of the proposed new quadrotor dynamic model identification method are promising.


INTRODUCTION
A wide area of robotics research is dedicated to aerial platforms.The quadrotor architecture has low dimensions, good manoeuvrability, simple mechanics and payload capability.
The study of kinematics and dynamics helps to understand the physics of the quadrotor and its behaviour.Together with modelling, the determination of the control algorithm structure is very important [1][2][3][4][5][6].The quadrotor unmanned aerial vehicle (UAV) is controlled by angular speeds of four motors.Each motor produces a thrust and a torque, whose combination generates the main trust, the yaw torque, the pitch torque, and he roll torque acting on the quadrotor.Motors produce a force proportional to the square of the angular speed and the angular acceleration; the acceleration term is commonly neglected as the speed transients are short thus exerting no significant effects.Soft computing methods can be efficiently applied together with and also instead of conventional controllers.
Fuzzy modelling [7][8][9][10][11][12] can be conducted as black-box modelling where all the system knowledge is mere input-output data, however when expert knowledge is readily available, we should take advantage of itfuzzy grey-box modelling is a rational choice.Identification of linear parameters is a well-studied area, with efficient matrix algebra and singular value decomposition based reliable tools.Non-linear parameters can also be simply traced to their local optimum with well-studied gradient descent methods, but we should always keep in mind that gradient descent methods are trapped by local optimum areas.Evolutionary algorithms are robust global optimum search engines, capable of multi-objective search as described in [13][14][15][16].
The article is organized as follows.In Section 1 the Introduction is given, in Section 2 the quadrotor dynamic model is presented.In Section 3 the Fuzzy-logic systems are illustrated.In Section 4 the multi-objective Genetic algorithms are illustrated.Section 5 presents the simulation setup and simulation results.Conclusions are given in Section 6.

QUADROTOR DYNAMIC MODEL
Motors of a quadrotor can only turn in a fixed direction, so the produced force is always positive.Motors are set up so that two opposites form a pair, which turns clockwise, while the other pair rotates counter-clockwise.This arrangement is chosen so that gyroscopic effects and aerodynamic torques are canceled in trimmed flight [17][18][19][20].
The main trust is the sum of individual trusts of each four motor.The pitch torque is a function of difference in forces produced on one pair of motors, while the roll torque is a function of difference in forces produced on other pair of motors.The yaw torque is sum off all four motor reaction torques due to shaft acceleration and blades drag.The motor torque is opposed by a general aerodynamically drag.
The complete dynamics of an aircraft, taking into account aero-elastic effects, flexibility of wings, internal dynamics of the engine, and the whole set of changing environmental variables is quite complex and somewhat unmanageable for the purpose of autonomous control engineering.
For a full dynamic model of a quadrotor system both (1) the center of mass position vector of (x, y, z) in fixed frame coordinates and (2) the orientation Euler angles: roll, pitch, yaw angles (Φ, θ, ψ) around body axes X, Y, Z are considered for the vector of generalized coordinates q.Using the Euler-Lagrange approach it can be shown how the translational forces F ξ applied to the rotorcraft due to main trust can be full decoupled from the yaw, pitch and roll moments For a full dynamic model of a quadrotor system both (1) the center of mass position vector of (x, y, z) in fixed frame coordinates and (2) the orientation Euler angles: roll, pitch, yaw angles (ϕ, θ, ψ) around body axes X, Y, Z are considered for the vector of generalized coordinates q.Using the Euler-Lagrange approach it can be shown how the translational forces F ξ applied to the rotorcraft due to main trust can be full decoupled from the yaw, pitch and roll moments τ: where m is the quadrotor mass and g is the gravitational constant.
where J is a 3  3 matrix, called the inertia matrix and C is also a 3  3 matrix that refers to Coriolis, gyroscopic and centrifugal terms.Further on, for the scope of this article we shall address only equation (2) as the quadrotor dynamic model to be identified.
Equation ( 2) can be analyzed as three resultant torques  i acting along the i th axes respectively as i ∈ (ϕ, θ, ψ), which using Christoffel symbols of the first kind can be defined as a function of the state vector of Euler angles q = (ϕ, θ, ψ), their velocities ( ) and accelerations ) as: The first component of equation ( 3) is shortly referred to as ̈ the inertia matrix part, while the second as C q  the Coriolis matrix term for which components are defined as: Where D ik , D ijk are in general, highly non-linear scalar functions of the state vector q.They contain sin( .) and cos( .) functions of q, and their products and sums defined by the geometry of the system.
There are general relations that can be used for reducing the number of unknown elements of J and C, like: (1) J is symmetric and (2) D ijk are Christoffel-symbols of D ij , thus further properties are inherently defined as: It should be noted that direct measurement of any single component from equation ( 4) is not possible; the only measurable data, on the output of the system, is the resultant torque of equation (3).
Identification of all non-linear functions (4) under these terms is a considerable problem.

FUZZY -LOGIC SYSTEMS
Takagi-Sugeno-Kang (TSK) type Fuzzy-logic systems (FLSs) having n inputs and 1 output are defined as: where M is the number of rules, q is the vector of n input variables, y l is a scalar function of n input variables, defined by (n + 1)c parameters as in equation ( 8).The antecedent, the premise part of a fuzzy rule is: is the membership function (MF) of the i th input variable in the l th rule that defines the linguistic value F l(i) .The linguistic form of the l th rule from the previously described first order TSK FLS is defined in [13] as: Zadeh-formed MFs are the z-, the s-, and  -functions (named after their shape) defined respectively as: where b 1  b 2  b 3  b 4 are parameters defining MFs.In case there is more than one value q such that the degree of membership of q is equal to one, the interval where A linguistic variable can be assigned K different linguistic values, each described by a MF , the MFs are said to form a fuzzy-partition.Forming fuzzy-partitions by antecedent membership functions ensures that there can not be a numerical input within the defined input range that will not result in firing at least one rule consequent of the fuzzy model, which means that there is a defined output for all possible input states.Keeping specific properties of fuzzy-partitions imposes a set of hard constraints on membership function parameters as detailed in [15].By imposing these restrictions on all linguistic variables of the FLS, and additionally assuming that the rule base is complete in the sense that it covers the whole input domain, it immediately follows that the TSK model structure of equation ( 6) simplifies to: Automatic fine tuning FLS parameters that satisfies all of above listed constraints is a significant problem.
In [15] a method is introduced that simplifies parameter optimisation of equation ( 11) while preserving all required constraints.Fuzzy-partitions can be simply formed from Zadeh-typed MFs by making equal the last two parameters of each preceding MF to the first two parameters of the succeeding MF.This way a fuzzy partition of K MFs is defined by 2(K -1) + 1 parameters.Let our input space be normalised (x min = 0 and x max = 1).If we do not want to allow any plateaux, parameter b 2 must be equal to b 3 in equation ( 9) this way the number of parameters is further reduced to K -2.When we take into consideration all constraints of equation ( 10), we end up with a series of strictly ordered parameters:  ,  ,  , ( , where b 0 = 0 and b K-1 = 1.This way the ordered series of K-2 parameters b i , together with constants 1 and 0, are the minimal number of parameters to define a fuzzy-partition of Zadeh-formed MFs.This minimal number of non-linear parameters is a very important issue for optimisation as over parameterised systems are hard to optimise.The only problem is that when tuning the non-linear parameters of a FLS having an n dimensional input space, we must comply with   pieces of hard constraints.Although there are a number of constrained optimisation methods it is obvious that an unconstrained optimisation method would be more efficient.Let us consider K -1 pieces of rational, positive or zero parameters as proposed in [12]: (17) When we define b k as: for every k = 1, …, K -2; all the constraints of equation ( 12) and equation ( 13) are automatically fulfilled for every a  from equation ( 18) without any further restrictions on any a  , other than 0  a  .An ANFIS like optimisation, defined in [16] or any other efficient unconstrained nonlinear numerical method can be applied to minimise equation (11) error along the  a parameters.For calculating all linear parameters a linear least square (LS) method can be applied to c l(j) parameters of the consequent part.To avoid traps of local optimal solutions for  a , a preliminary global search should be applied.

MULTI-OBJECTIVE GENETIC ALGORITHMS
A genetic algorithm (GA) is constructed on bases of imitating natural biological processes and Darwinian evolution [21][22][23][24].GAs are widely used as powerful global search and optimization tools [25].Real life optimization problems often have multiple objectives.To establish ranking of chromosomes for Gas the comparison of two objective vectors is required.Often a simple weighted sum is used, but its drawbacks are widely known.Pareto based comparison [19] is the bases of a few popular methods like Non-dominated Sorting GA (NSGA) [22] and Multi-Objective GA (MOGA) [23,25].A general multi-objective optimization problem consists of n number of scalar minimization objectives where every scalar objective function f i (x) is to be minimized simultaneously, where x is an n-dimensional vector of parameters.As maximization can be easily transformed to minimization, the generality of the previous statement stands.A vector x 1 Pareto-dominates x 2 , when no scalar component of x 2 is less than the appropriate component of x 1 , and at least one component of x 1 is strictly smaller than the appropriate component of x 2 .Since no metrics can be assigned to Pareto-dominance, there have been two different attempts to define a GA ranking method, which can be used for Pareto-dominance vector comparison: (1) "Block-type" ranking is defined in [23] as: Rank is equal to 1 + (number of individuals that dominate the i th individual); (2) "Slice-type" ranking is defined in [5] as: Rank is equal to 1 + (number of turns when the non-dominated individuals are eliminated, needed for the i th individual to become non-dominated).
Quantity-dominance, as proposed in [15] is defined as: vector a = [a i ].Quantity-dominates vector b = [b i ] if a has more such a i components, which are better than the corresponding b i component of vector b, and a has less such a j components, which are worse than the corresponding b j .A metrics can be defines as: the measurement of the extent of Quantity-dominance is the difference between the number of better and the number of worse components.For a measurement based ranking method the Rank of the i th objective vector can be simply defined as the sum of Quantity-dominance measurements for every individual measured from the i th vector.This ranking method can be readily applied with Quantitycomparison.The proposed vector comparison method provides more information when comparing two vectors than the classic Pareto-based comparison, thus the GA is faster, more efficient in its search.The MMNGA algorithm is computationally less expensive, and more efficient compared to the classical methods, and its results analyzed on a number of GA hard problems are at least equally good [16].
In case of multi-rotors the roll and pitch are equal to: . sin cos atan , ) ( cos sin asin From equations ( 3) and (19) it is obvious that controll torques for multi rotors are direct functions of up to the forth time derivatives of state variables (x,y,z) and ψ.To have realistic, feasible torques along a trajectory, which are efficiently controllable without chattering, we need smooth torque changes.Having  = (, ̇, ̈) and  = (, ̇, ̈) for smooth torque changes, we need smooth changes of the so called displacement crackle () =  5 / 5 , the fifth time derivative of displacement r(x, y, z).Proposal of this article is to use a smooth displacement crackle function, which can be defined with a continuous displacement pop function () =  6 / 6 as: (20) where  d is the natural dampened frequency from equation (3),  d the phase delay is kept is zero, and the gain G is selected for each trajectory and system such that the required boundaries for displacement, velocity and acceleration limits are met.The integration constant for () is to be selected as equal to  to achieve the required properties for the cracle function; for all further integrations to calculate trajectory snap, jerk, acceleration, velocity and displacement by intergrating () we are to use integration constants equal to zero.The resulting general trajectory plot is as presented in Figure 1.
We can eficiently identify D ij inertia matrix components of the dynamic model in equation ( 4) as FLSs defined by equations (11) to (18), where the FLS general input variable q will be substituted for the appropriate state variables of (ϕ, θ, ψ).When the D ij inertia matrix components are constructed in this way, forming the D ijk components as Christoffel symbols is to be expressed by partial derivatives of equation (11) : The unknown inertia matrix components of equation ( 2) to be identified are: ) , ( ) Based on quadrotor system structure and inertia matrix symmetry the remaining inertia components are known to be: Remaining D ijk components are trivial identities defined by equation ( 5).

SIMULATION SETUP AND RESULTS
The proposed method is tested for a quadrotor system simulation from [1] with following parameters: gravity constant g = 9,81 m/s 2 , mass m = 6 kg, trust factor k = 121,5 e -6 , drag factor b = 2,7e -6 , body inertia along axes X, I X = 0,6 kgm 2 , body inertia along axes Y, I Y = 0,6 kgm 2 , body inertia along axes Z, I z = 0,6 kgm 2 , simulation time t = 3 s.The training data set is collected from a simulation along a trajectory with jounce for (x,y,z) and ψ defined so that position changes simultaneously along a unit cube main diagonal (0 0, 0)-(1, 1, 1), while performing full circle rotation in jaw motion 0-2.The simulated resultant torque training data is as presented in Figure 2. The calculated roll, pitch and yaw motions are presented in Figs.3-5.
Non-linear a K parameters of the system are identified in a manner that first the input space is normalised to the unit hyper-cube.A set of non-linear parameters consists of six times four a K integer parameters for defining six fuzzy-partitions of five MFs each, where each partition consists of one z-type MF, three π-type MFs and one s-type MF as in ( 9)- (18).These six fuzzy-partitions serve as antecedents for the four fuzzy-systems like in equation ( 11) and ( 21), used for identifying D ij , with ij = (13,22,23,33) as defined in equations ( 22)-( 24) and (5).Two unknown linear parameters D 11 and D 12 of the quadrotor model as in equation ( 23), together with 170 linear parameters of the four TSK FLSs (2 FLSs with 5 MFs on one input, each rule with 2 c parameters, plus 2 FLSs with 5 MFs on both of the 2 inputs, each rule with 3c parameters) of equations (22) and equations (24) are determined by the SVD-based LS method as used in [15].Concluded from equation ( 17) six fuzzy-partitions (antecedent part of 2 FLSs with 1 input, plus 2 FLSs with 2 inputs are covered by 6 independent fuzzy-partitions) are represented by a vector of six times four K a parameters, which are optimized by a multi-objective hybrid genetic algorithm as detailed in [16].Each chromosome evaluation is extended to include an additional round of nonlinear LSQ optimization of K a parameters.Chromosomes are updated before applying further GA operators, so the GA does not waste time on local optimization; only global search capabilities of the GA are utilized.Three objectives are set for minimization: (1) the root mean square of the torque identification error, (2) the maximum absolute error for any given training data input-output pair, and (3) the condition number of the linear system of equations used for LS calculation of linear parameters.The GA is set to work on a population of 125, divided into 5 subpopulations with migration rate 0,2 taking place after each 5 completed generations.Crossover rate, generation gap and insertion rate is 0,8, selection pressure is 1,5.In each generation 4 % of individuals are subject to mutation, when 1 % of the binary genotype is mutated.Individuals, chromosomes are comprised of 24 Gray-coded integers, each consist of 16 bits.The initial population is set up in a completely random manner.Matrix of the linear equation is preprocessed from equation (3), where FLSs like equation (11) and their partial derivatives like equation ( 21) are inserted as described in equations ( 22)- (24).Unknown linear parameters are D 11 , D 12 and the 170 c parameters of fuzzy-rule consequents.
Evaluation of each individual is done as follows: (1) Convert the coded a K values from the chromosome to b k by equation ( 18).(2) Evaluate all MFs, which will comprise six fuzzy-partitions from each of six b k quadruplets by equations ( 14)- (16).Also evaluate derivatives of equations ( 14)-( 16).(3) The pre-processed matrix of the linear equation is evaluated with the MFs.(4) Linear components of equations ( 11) and ( 21) are calculated by SVD decomposition as described in [20].(5) Next the K a parameters are fine-tuned by the Matlab "lsqnonlin" function, (6) MFs are re-calculated for the optimised a K parameters and all linear parameters are re-calculated.(7) Resulting optimised K a parameters are re-assigned into the chromosome of the evaluated chromosome.For the multi-objective rank assignment described in [16], the objective vector is created from: (i) the mean square of the identified torque error, (ii) the maximum absolute torque identification error and (iii) the condition number of the matrix of the linear equation.Stochastic universal sampling is used for selecting the next generation without explicit elitism.To speed up the GA processing, a database of evaluated chromosomes and their objective vectors is created, so only unique new individuals are evaluated in each generation.Convergence is achieved in some 50 generation evaluations, when the mean square error is in order of 5e -7 , the maximum torque error is smaller than 0,005 Nm.For non-dominated chromosomes the condition number of the used matrix of linear equation is in order of 1e

CONCLUSIONS
Simulation results of the proposed new quadrotor dynamic model identification method are promising.The quality of identification with the relative torque error being uniformly <0.1% is excellent, suitable for taking part in a model based control algorithm.The typical condition number for used linear parameter evaluations is very high for the used training data setup, so a more advanced trajectory has to be planned with sufficient inertia excitation along the complete input domain.Also the FLS structure is to be made flexible, in terms that the GA should be able to turn off unnecessary MFs and thus reduce the number of unnecessary rules and linear parameters.
called the plateau of the  k MF.When having for example 3 naturally ordered linguistic values l  {a, b, c} (a  low, b  medium, c  large) constraints on parameters to preserve this ordering are: ) Let us add two more constraints: 0 < b 1 and b K-1 < 1. (13) Let us define the first MF to be: the K-th, the last MF concluding the fuzzy partition be: Let us define intermediate kth MFs to be:) , on equation (5) the following Coriolis term matrix components can be calculated by equations (

Figure 3 .
Figure 3. Roll training data for input.

Figure 4 .
Figure 4. Pitch training data for input.

Figure 5 .
Figure 5. Yaw training data for input.