Combination of terminal sliding mode and finite-time state-dependent Riccati equation: Flapping-wing flying robot control

A novel terminal sliding mode control is introduced to control a class of nonlinear uncertain systems in finite time. Having command on the definition of the final time as an input control parameter is the goal of this work. Terminal sliding mode control is naturally a finite-time controller though the time cannot be set as input, and the convergence time is not exactly known to the user before execution of the control loop. The sliding surface of the introduced controller is equipped with a finite-time gain that finishes the control task in the desired predefined time. The gain is found by partitioning the state-dependent differential Riccati equation gain, then arranging the sub-blocks in a symmetric positive-definite structure. The state-dependent differential Riccati equation is a nonlinear optimal controller with a final boundary condition that penalizes the states at the final time. This guides the states to the desired condition by imposing extra force on the input control law. Here the gain is removed from standard state-dependent differential Riccati equation control law (partitioned and made symmetric positive-definite) and inserted into the nonlinear sliding surface to present a novel finite-time terminal sliding mode control. The stability of the proposed terminal sliding mode control is guaranteed by the definition of the adaptive gain of terminal sliding mode control, which is limited by the Lyapunov stability condition. The proposed approach was validated and compared with state-dependent differential Riccati equation and conventional terminal sliding mode control as independent controllers, applied on a van der Pol oscillator. The capability of the proposed approach of controlling complex systems was checked by simulating a flapping-wing flying robot. The flapping-wing flying robot possesses a highly nonlinear model with uncertainty and disturbance caused by flapping. The flight assumptions also limit the input law significantly. The proposed terminal sliding mode control successfully controlled the illustrative example and flapping-wing flying robot model and has been compared with state-dependent differential Riccati equation and conventional terminal sliding mode control.


Introduction
The sliding mode control (SMC) and variable structure systems are among the most popular nonlinear robust methods; they find the control law by defining a proper sliding surface. The SMC has so many branches such as a second-order, adaptive, terminal, and also in combination with other techniques including fuzzy, neuralnetwork, optimal, optimal-integral, etc. 1 The focus of this current research is on terminal sliding mode control (TSMC). The TSMC is an updated structure of the SMC, introduced in the 90s. 2 The terminal attractors in the new controller guarantee the finite-time state convergence. The nonlinear definition of the sliding surface makes the manifold an attractor. 3 Yu et al. 4 used a continuous TSMC for controlling rigid robotic manipulators. The TSMC was improved in many aspects, such as a nonsingular version of the controller, 5,6 fast TSMC, 7 fractional-order one, 8 adaptive nonsingular, 9 etc. Wang and Xu 10 investigated the micro-positioning control system via an adaptive TSMC. The micropositioning system used piezoelectric actuators that had resonance; the adaptive novel gain controlled the uncertain system with high performance. Xiong and Zhang 11 employed global fast TSMC for quadrotor control in set-point regulation. The system was divided into two fully-and under-actuated system and controlled successfully.
This current work intends to find a finite-time controller with robust characteristics. The TSMC is a good selection though the final time cannot be set as a control input parameter. 12 Shiri et al. 12 solved the problem by dividing the sliding surface into two linear (conventional SMC) and nonlinear sections (nonsingular TSMC). Then it was possible to define the final time as an input to the control system. The controller was applied for a van der Pol oscillator and a robotic manipulator. Here in this current work, the time setting approach is used and the final time could be defined as an input to the control system. On the contrary to Shiri et al., 12 here in this work, one nonlinear sliding surface is considered. To have a command over the final time, an optimal gain of the state-dependent differential Riccati equation (SDDRE) is used. The SDDRE is a nonlinear closed-loop suboptimal controller with the power of finite-time control by penalizing a terminal boundary condition. 13,14 Tripathy et al. 15 used the SDDRE for manipulator control. There are several solution methods to the differential Riccati equation, such as backward integration, 16 Lyapunov-based method, 17,18 and state transition matrix approach. 19 Here we used the Lyapunov-based approach as the solution to the SDDRE with negative gain. The gain of the SDDRE could be partitioned and used in other controllers preserving the finite-time characteristics.
Recently, the flapping-wing flying robot (FWFR) has been a challenging model for the assessment of the controllers in terms of robustness and performance. The model is highly nonlinear and coupled with force/ moment components. One of the open problems in this field is the dynamic modeling of a flapping-wing system. The six-degree-of-freedom (DoF) in modeling leads to complicated systems especially in terms of lift, drag, and aerodynamics calculations. The modeling of aircraft 20,21 is the base of the modeling though the flapping wing adds complexity and changes to the system. The dynamic derivation of the gliding phase of the flapping systems is very similar to aircraft dynamic; however, the inertia and weight, and propulsion forces are significantly lower. The investigation began by looking into the bird's flight dynamics. 22, 23 Taylor and Thomas 22 studied the flapping flight dynamics and longitudinal stability of animal flights. Deng et al. 24 researched flapping flight dynamics for insect-sized systems and studied the lift, drag, and body dynamic. Platzer et al. 25 reviewed the progress and challenges of flapping-wing flight considering different wing configurations. It was found that thrust, lift, and propulsive efficiency on flapping mode depend on the frequency and amplitude of the flapping. The insect-sized flapping-wing systems have been attractive since the actuation of them and the weight of the system provided the flight situation more comfortable than large-scale flapping systems. [26][27][28][29][30][31][32][33] In that case, significant load-carrying capacity is not expected.
Increasing the size of the FWFR requires a complicated mechanical mechanism to obtain constant flapping and control over the frequency. Armanini et al. 34 used flight data of experiments to define ornithopter dynamics through global system identification. Hassan and Taha 35 showed that the flight path of the flapping robot was subjected to vibration. This result was found by relaxing the common assumption that neglected the wing inertial effects and directly averaged the dynamics over the flapping cycle. Bakhtiari et al. 36 studied the modeling focusing on the aerodynamics effects on the edge of the tail, experimentally. All the dynamics modeling (some works with the experimental investigation on lift, drag, and aerodynamics) and control implementations (mostly in theoretical schemes and simulations) lead the researchers toward carrying more weight to apply the FWFR for a particular task. [37][38][39][40] One application is perception and visual inspection by FWFR which requires more lift capacity due to necessary external devices, such as a camera, digital board, and heavier batteries. 39 The majority of the modeling used general dynamic modeling of a rigid body in gliding conditions with added lift, drag, and aerodynamics forces. The representations of the rigid body were similar though the aerodynamics modeling was different based on the situation and condition. In this work, six-DoF modeling of the FWFR using rigid body dynamics is considered. The aerodynamics of the system is added to the dynamic equation and the relation between the force/ moment inputs and rudders and thrust actuators is presented using Taylor series expansion (linearized model). It should be emphasized that the model is nonlinear and coupled, and the relation between force/moment and actuators is only linearized.
The main contributions are as follows: C1. Introducing a terminal sliding mode to control a class of nonlinear uncertain systems in a predefined finite-time; that is defined as an input parameter; C2. Stability proof of the TSMC gain for a predefined finite-time control using the partitioned SDDRE control gain; C3. Six-DoF fully actuated SDDRE-based TSMC control of FWFR (modeling and simulation).
The rest of the article is as follows: Section ''TSMC'' expresses the control structure of the proposed TSMC. Dynamic modeling of the FWFR is expressed in section ''Dynamic modeling of FWFR.'' State-space form and control implementation are reported in section ''State-space form and implementation.'' Section ''Simulations'' presents simulations: illustrative example, FWFR, and comparison with conventional SDDRE and TSMC. Concluding remarks are reported in section ''Conclusion.'' Notations R n denotes the n-dimensional Euclidean space; R n 3 m is the set of m 3 n real matrices; (Á) T is the transpose of a matrix or a vector; diag(Á) means a diagonal matrix; I n 3 n and 0 n 3 n denote n 3 n identity and zero matrices. j Á j defines absolute value and jj Á jj 2 is norm-2 of a vector or a matrix. eig½A provides a diagonal matrix with eigenvalues of A.

System definition
Consider a second-order nonlinear differential equation (1) presents a fully actuated mechanical system. The generalized coordinates of the system (equation (1)) are collected in q 2 R n vector and the input force or torque in u 2 R n in which n denotes the degree-of-freedom (DoF). The system is limited to a common form that covers a wide range of mechanical systems by defining f(q, _ q) =À M À1 (q)(c(q, _ q) + g(q) + D(q, _ q) _ q) and E(q) = M À1 (q) where M(q): R n ! R n 3 n is inertia matrix, c(q, _ q): R n 3 R n ! R n is Coriolis and centrifugal vector, g(q): R n ! R n is a vector that includes gravity terms, and D(q, _ q): R n 3 R n ! R n 3 n collects all terms related to friction, drag (underwater or flying systems), aerodynamics forces, etc.
Selecting the state vector as x 2n 3 1 = ½q T , _ q T T , the reduced-order differential equation (state-space representation of equation (1) where a(x) : R 2n ! R 2n and g(x, u): R 2n 3 R n ! R 2n are piecewise-continuous vector-valued functions that satisfy the Lipschitz condition.
The symmetric gain of state-dependent Riccati equation In this section, the finite-time suboptimal gain of the state-dependent Riccati equation (SDRE) will be partitioned into four square blocks and will be used in the definition of the sliding surface in section ''The SDRE-TSMC design.'' The SDRE control law will not be used directly; however, the finite-time gain will be used in the definition of the TSMC sliding surface. The first step is to rewrite state-space representation (equation (2)) in the state-dependent coefficient (SDC) parameterization 0 n 3 n j I n 3 n À À À j À À À À À À À À À À À À À Remark 1. The gravity vector usually contains trigonometric functions (robotic arms) or constant values (flying objects). The cos (Á) function cannot be factored after using Taylor series expansion since the first term is 1. Also dividing cos (Á)=x i imposes singularities at equilibrium x i = 0. So, the gravity vector will be removed from SDC matrices (equation (3)) and will be compensated in TSMC control law.
Condition 1. fA(x), B(x)g is a completely controllable parameterization of the system (equation (2)) for all x in t 2 ½0, t f . 41  The cost function is where R(x): R 2n ! R n 3 n is input weighting matrix, symmetric positive-definite, and Q(x): R 2n ! R 2n 3 2n is state weighting matrix, symmetric positive-semi-definite. F 2 R 2n 3 2n also penalizes the states at the terminal time t f ; F = F T and F ø 0.
Condition 2. fA(x), Q 1 2 (x)g, in equation (4), is a completely observable parameterization of system (equation (2)) for all x in t 2 ½0, t f , where Q 1=2 (x) is Cholesky decomposition of Q(x). 41 Defining all matrices, A(x), B(x), Q(x), R(x), and F, one could solve the SDDRE with terminal boundary condition K(x(t f )) = F. Lemma 1. The nonlinear system (equation (3)), with SDC pair based on Conditions 1 and 2, can be stabilized using the standard SDDRE control law is symmetric positive-definite solution to SDDRE (equation (5)). 42 The proof of Lemma 1 could be found using the Lyapunov candidate V(x) = x T K(x)x. The details of solution methods to equation (5) could be visited in Korayem and Nekoo;42 we suppose that SDDRE could be solved with one of the methods such as the Lyapunov-based approach with negative gain.
The control gain, the solution to the SDDRE, is partitioned into four square blocks 43 and is substituted in overall control gain in which The gain X 1 (x) is multiplied by the position vector, q(t), and X 2 (x) by the velocity vector, _ q(t). From a general point of view X 1 (x) acts as proportional gain and X 2 (x) as a derivative one.
Remark 3. The proportional and derivative gains X 1 (x) and X 2 (x) possess finite-time characteristics for the regulation of the error and velocity of the error vector to zero, tuned by F.
The finite-time proportional gain of the SDDRE is required for the definition of the sliding surface in section ''The SDRE-TSMC design,'' though it is not a symmetric positive-definite gain. To define the gain in that shape and preserve the finite-time characteristics, it will be transformed into 43 where K SP (x) is made symmetric by multiplication of K 12 (x)M À1 (x) from the left-hand side of the equation and leveled by division to the norm-2 of that K 12 (x)M À1 (x) 2 . If the weighting matrix of states (preserving symmetric positive-definite condition) at final time is partitioned into four submatrices then, the F 12 sub-block will define the final boundary condition of K 12 (x) in equation (6), as K 12 (t f ) = F 12 . However, F 12 is sign-indefinite and for bounding K SP (x) in section ''The SDRE-TSMC design,'' a positivedefinite F 12 is required. As a result, the following conditions are considered: F 11 , F 12 , and F 22 are diagonal matrices with positive components, and F 12 \ ffiffiffiffiffiffiffiffiffiffiffiffiffi . This provides a positive-definite symmetric F and also F 12 .
Assumption 1. It is assumed that F 11 , F 12 , and F 22 , in equation (7), are diagonal matrices with positive components, and F 12 \ ffiffiffiffiffiffiffiffiffiffiffiffiffi Assumption 1 implies that the final weighting matrix must be chosen relatively bigger than the weighting matrix of states, F ) Q(x), or in other words, F 12 ) ½Q 12 (x) n 3 n to force the solution to start from an initial condition and move to a large final boundary condition.

The SDRE-TSMC design
Consider the state-space representation of the general nonlinear affine-in-control time-invariant system (equation (2)). The goal is to define a robust nonlinear control law to regulate the error of the system to the desired condition in predefined terminal time. The terminal sliding mode is a finite-time controller in nature though the user cannot set the final time as an input to the control problem using common standard forms. The terminal sliding mode controller is finite time in nature; however, the final time cannot be introduced as an input to the control system. The reason is that the total final time is the summation of rise time and settling time, t f = t r + t s . Calculation of t s needs the computation of e(t r ), and this value (error at rise time) is unknown before the simulation of the system or implementation of the experiment. So, while the terminal time is finite, the user does not know the exact final time before the implementation. Shiri et al. 12 used two sliding surfaces to present a structure accepting the final time as an input to a nonsingular terminal sliding mode controller. A linear sliding surface was responsible to define the reaching time and the nonlinear one set the convergence time; together they could determine the total final time. The combination of terminal sliding mode and finite-time state-dependent Riccati equation provides a solution for determining the final time before the implementation of the controller and forces the system to finish its control task by weighting matrices of the SDRE. Here, one nonlinear sliding surface is selected with a finite-time gain to present a simpler structure.
The nonlinear sliding surface is designed as for i = 1, . . . , n; b . 0, q and p are positive odd integers following p . q.
Remark 4. The definition of the sliding surface uses the relative part to position variables of the system in regulation, K SP (x). Observing the selection weighting matrices for the SDRE in the literature, 42,[44][45][46] it can be seen that for precise and fast position control using optimal control, the velocity states must be relaxed. So, the role of the proportional gain is more highlighted. To keep the focus on position states and also present a sliding surface as simple as possible, the control gain related to derivative states, X 2 (x), is not included in equation (8).
The error vector is defined as e(t) = q(t) À q des (t) in which q des (t) 2 R n is the desired values for the generalized coordinates. Derivation of the error also results in _ e(t) = _ q(t) À _ q des (t) in trajectory tracking and _ e(t) = _ q(t) in regulation; and Assumption 2. The unknown system is bounded above, the vector f i (x) łf i (x) and matrix P n Assumption 3. The deviation of the external disturbance, f(x) and E(x), from their bounds is limited to the known values The convergence and stability of the TSMC have two parts: Stage 1, showing the convergence of the system to the sliding surface, and Stage 2, the convergence of the error to zero on the sling surface. 5,47 The latter is done by computation of the settling time on the sliding surface s = 0. The both steps could be expressed in one theorem, such as Feng et al. 5 and Xu et al.; 47 however, for simplification in finite-time calculation, we presented it in two Theorems 1 and 2. The stability discussion of the proposed TSMC is studied by the following steps: First, the convergence of the system to the sliding surface in reaching time, t r , will be shown and proved, Theorem 1. Then the computation of the reaching time is presented in Lemma 2. Finally, the convergence of the error (along with the sliding surface s = 0) to zero in finite time will be shown by computation of the final time, T = t r + t s , where t s is the settling time, Theorem 2.
Theorem 1. (Convergence to the sliding surface). Considering Assumptions 2 and 3, the nonlinear system (equation (2)) will reach the sliding surface (equation (8) and e s are set along with the proof.
Proof. A Lyapunov candidate is used V = 0:5 s T s to guarantee stability through its derivative _ V = s T _ s \ À h T jsj, known as the sliding condition. The derivative of the sliding surface (equation (8)) should be computed (please see the details in Appendix 1) to form the condition _ s i = 0, for providing the equivalent control part. Substituting € e i , equation (9), in equation (13) and mathematical manipulation, one could present is replaced with the upper bounds of the system dynamics, Assumption 2, Expressing the matrix/vector presentation of equation (14), the equivalent or nominal control is presented where e s := . . .
A correction part will be added to the nominal control (equation (15)) to complete the TSMC input law The gain K T (x): R 2n ! R n 3 n is determined by the sliding condition s T _ s \ À h T jsj. Substituting the actual system (equation (1)) and control law (equation (17)) into sliding condition, it is rewritten as The sliding condition (equation (18)) is represented in compact form The scalar form of equation (19) is which can be represented as Dividing equation (20) by s i removes i-summation then computing the absolute value of that results in which could define Changing equation (21) to matrix form will obtain a vector Substituting G(x) in equation (22) and considering Assumption 3 result in inequality in terms of the defined bounds, using equations (10)- (12) K represents the diagonal arrays of K T (x) = diag( K T (x)). Note that equation (12) results in jE min (x)Ê(x)j À1 ø jE(x)j À1 . j Time analysis is a critical subject in TSMC, and computation of the settling time is done in the literature. The required time for the error, e(t), to move between the initial condition and the reaching point of the sliding surface is the so-called reaching time, t r ; and the remaining time for the error to converge from reaching time e(t r ) to equilibrium point (with a specified error band) is the so-called settling time t s . The sum of reaching and settling time provides the necessary time for convergence of the TSMC T = t r + t s (terminal time), see Figure 1. In this work, the necessary time of convergence, T, is different from final time of operation, t f . Here we define t f as input parameters; however, T might be bigger or less than that, and is known after the simulation or experiment based on the error value. Usually, the control designer prefers a faster final time.
Lemma 2. (Reaching time) 48 The reaching time of ith error e i (t) to sliding surface is defined by t r, i ł js i (0)j=h i .
Regarding the sliding condition in Theorem 1, the reaching time of ith error in Lemma 2 is computed by integrating sliding condition 0:5d=dts 2 i \ Àh i js i j between 0 and t r, i : Lemma 3. (Settling time of conventional TSMC) 49 The settling time of the conventional TSMC is defined by t s, i = je i (t r, i )j 1À q p =b(1 À (q=p)). At reaching time, t r, i , the sliding surface is zero which results in de i =dt =À be i je i j q p À1 , and after integration results in t s, i = je i (t r, i )j 1À q p =b(1 À (q=p)) for conventional TSMC.
If the gain K SP (x) is removed from sliding surface (equation (8)), the control structure will be similar to a conventional TSMC. So, the gain K SP (x) manipulates the time in control structure. To compare the final time of the control problem with conventional TSMC, Lemma 3 has been expressed (in detail) for comparison.
Theorem 2. Considering Assumption 1, eig½K SP (x(t)) ł F 12 , the error of system (equation (2) Proof. The reaching time is defined based on Lemma 2 and the value of settling time t s, i will define the terminal time T. At reaching time, t r, i , the sliding surface is zero We consider the matrix-vector form of equation (23) for the sake of simplicity and rewrite it as The upper bound of K SP (x) will be considered to compute the integral where K C = F 12 is the constant upper bound of eig½K SP (x) ł K C , stated in Assumption 1. Considering K SP (x), the integration of equation (24) is not possible since the gain is a time-varying state-dependent value (related to error variable). We consider a constant gain K C for integration, then equation (25) is integrated as Considering equation (26) and if time tends to the final time of control task, t f , the control gain goes to final boundary value, lim t!t f K SP (x(t)) = F 12 . Each row of right-hand side of equation (26) which defines the terminal time T = max (t r, i + t s, i ) and that concludes the proof. It perfectly makes sense; if one increases F 12 , it will force the system to a faster response, see equation (27), and results in a shorter t s, i .
Remark 6. Theorem 2 expresses that the error of the system will converge to zero in finite time T, and this value is the shortest possible time, best case scenario. Considering putes a more realistic value, longer though finite, for the final time of error convergence.
Remark 7. The role of control over final time is played by the SDDRE control gain K SP (x) in equation (6), tuned by weighting matrices Q(x), R(x), and F. Remark 8. The role of chattering and elimination of that will be discussed in the ''Simulations'' section. To keep the derivation intact based on the sign(Á) function, saturation function was not considered in the derivation. However, to remove the chattering, saturation function or other similar ones, such as tanh ((Á)=e), could be used instead of sign(Á) in implementation.

Dynamic modeling of FWFR
The weight of the body is light and the contribution of the aerodynamics of the wing is significant to the model. In reality, birds can perform various stable maneuvers during flight, such as hovering, backward flight, diving, surging, etc. Such complete control over the motion expresses hybrid complex dynamics and continuous non-affine actuators. The mechanical design of a bird with all those properties might not be possible at the moment. This section intends to present a simplified six-DoF model of an FWFR with a minimum requirement of fully actuation mode. The coordinates and schematic model of the flapping-wing robot are presented in Figure 3. The FWFR is not limited to  a fully bioinspired design; so, a rudder was set on the tail for better control over the roll and yaw.
To cover a fully actuated model, six inputs are necessary, defined in Figure 3. Why six inputs? There are several flapping robots controlled experimentally by only two inputs on their tail and one input thrust (generated by constant flapping). 37,38,40 These models presented underactuated systems subjected to many limitations with a lack of control over the motion. The landing was also like a fall by reducing the flapping frequency. Considering six actuators simplifies the control problem and extends the model to be used by so many methods, and closes the gap between the mechanical models and real birds. The three additional inputs are folding parts near the tip of the wings represented by d RW and d LW for the right and left wing; and a unified rudder on the back of the wings to enhance the lift. Observing the flapping birds, it can be realized that they are moving continuously possessing curved shapes and different parts for changing the lift force on the wing. So, the selected inputs are u(t) = ½d RW , d LW , d L , F T , d UD , d LR T where d (Á) (rad) represents the rudders and folding blades on the tails and wings and F T (N) is the thrust force caused by flapping of the wings.
The following assumptions are held for the modeling and behavior of the system: Specifically, x c (t f ) should be significantly greater than y c (t f ) and z c (t f ). 4. The FWFR motion avoids severe changes in orientation and the main changes are limited to small deviations due to the flapping.
The fixed reference frame is set by (X, Y, Z), known as Earth frame; and the body (moving) coordinate is also presented by (x c , y c , z c ) placed on the center-ofmass (CoM) of the system. Four main forces are acting on the robot, the total weight of the system W, lift, drag, and thrust force. The generalized coordinates of the system are translation and orientation of the system in six-DoF, presented by q(t) = fj 1 (t), j 2 (t)g where j 1 (t) = ½x c (t), y c (t), z c (t) T (m) and j 2 (t) = ½f(t), u(t), c(t) T (rad). On the body frame, local linear and angular velocities are named as fy 1 (t), y 2 (t)g where y 1 (t) = ½u(t), v(t), w(t) T (m=s) and y 2 (t) = ½p(t), q(t), r(t) T (rad=s). The kinematics relation between the global and local frame, considering Flight Assumption 4, is presented 50 where in which, for example, c u = cos(u(t)) and t u = tan(u(t)).
The matrix T(j 2 ) = W À1 (j 2 ) holds in which and W(j 2 ) is invertible if det (W(j 2 )) = cos u 6 ¼ 0. The kinematics, equations (28) and (29), and the reference frames are in the common forms of flying objects in the aerospace community, such as aircraft, quadcopters, etc. The dynamics of the FWFR is where M(q): R 6 ! R 6 3 6 is the inertia matrix, ½C(q, _ q) _ q : R 6 3 R 6 ! R 6 is Coriolis and centrifugal forces, ½D(q, _ q) _ q : R 6 3 R 6 ! R 6 collects drag and aerodynamics terms, g(q): R 6 ! R 6 is gravity vector, and U 2 R 6 is input vector: in which m (kg) is the total mass of the FWFR, g (m=s 2 ) is gravity constant, and J j 2 ð Þ= W T j 2 ð ÞIW j 2 ð Þ, The drag matrix in the translation part is defined as where A i is the effective area of the system in i axis and C Di is the drag constant for i = x, y, z. Based on Flight-Assumption 4, the drag and aerodynamics effect in orientation dynamics are neglected.
The input vector includes two force and torque parts: The force/torque vector U has a nonlinear complicated relation with ultimate input vector u, which is nonlinear non-affine and unknown Representation of the equation (31) using Taylor series expansion around zero is All the rows of vector F(q, _ q, u) are multiplied by u, or in other words, F(q, _ q, 0) = 0. Considering the equilibrium F(q, _ q, 0) = 0 and neglecting the higher terms of the Taylor series, equation (32) is approximated by  where, that is, X d RW = ∂X=∂d RW : We call the coefficient matrix M x = ∂F(q, _ q, 0)=∂u a mixer matrix. In reality, for example, considering a real bird, it could be expected that all the inputs u contribute to the force/ torque vector U and all the elements of the mixer have values; however, this model is simplified by several assumptions. To find the effective elements of the mixer, the effect of the inputs is considered in decoupled mode, then based on Figure 3, one could rewrite equation (33) as The oscillation effect of flapping is more visible in z c axis and pitch angle u. Based on the flapping oscillation, a disturbance vector is applied to the dynamics, added to the system in state-space representation. The inputs are also bounded between d min ł d (Á) ł d max and 0 ł F T ł F T, max .

State-space form and implementation
The state variable of the system is selected as Considering the state vector (equation (34)) and the dynamics (equation (30)), the state-space form is obtained where D F = d F 3 ½0, 0, cos vt, 0, 0, 0 T is the disturbance caused by flapping in which d F is the amplitude of the oscillatory disturbance. Comparing equation (35) with equation (2), the val- ) and E(q) = M À1 (q) are defined. The SDC parameterization matrices are also shaped as À À À j À À À À À À ÀÀ À À À j À À À À À À ÀÀ which results in the gain K SP (x) = diag(K SP t (x), K SP o (x)). The nature of the orientation dynamic is faster than the translation. To solve the Riccati easier and avoid setting the eigenvalues of the dynamics near the imaginary axis, two separate SDDRE control is used.
It should be noted that considering the smooth and small variation of orientation variables, R ZYX (j 2 ) and T(j 2 ) possess almost identity matrices. That is necessary to guarantee the controllability and observability conditions, Conditions 1 and 2. The rest of the design is based on section ''The SDRE-TSMC design.''

Illustrative example
A van der Pol oscillator is simulated and compared with two other controllers. The state-space representation of the disturbed van der Pol equation is: 12 where C(t) = rand(1) 3 0:2 sin 100t is an external disturbance in which rand(1) generates a random number between 0 and 1 at each time step of the simulation. The proposed method is finite time and robust. It is good to compare with the SDDRE controller which is finite time though is not robust and also with the conventional TSMC which lacks optimality. The initial condition is x(0) = ½2, 1 T and the system is to be regulated to the equilibrium point. The control parameters are selected as R = 1, Q = I 2 3 2 , F = 10 3 Q, h = 2, b = 1:5, q = 7, and p = 9.
We also considerf(x) = 1:2 3 (À x 1 + (1 À x 2 1 )x 2 ), E = 1:1, and C c (t) = 0:2 sin 100t which defines (e.g. here we considered rand(1) = 0:501) the constants, based on equations (10)-(12): The parameters of the SDDRE and TSMC controllers are similar weighting matrices and parameters. The only change in the TSMC is the constant gain of K SP = 1. The comparisons have been done in different terminal times to illustrate the role of predefined time setting in the results. The error and input of the system in different final times are presented in Figure 4, and the error details are in Table 1. The error of the proposed TSMC was obtained less than the TSMC and SDDRE which showed for a system with uncertainty, a combination of the robust and optimal control could score better results. For the 6s final time, the proposed TSMC improved (Possible improvement (%): ((old data Ànew data)=old data) 3 100) w.r.t. the TSMC 61% and SDDRE 61:3%; for the 4 s final time, the proposed TSMC improved w.r.t. the TSMC 56:3% and SDDRE 87:3%; and for the 2 s final time, the proposed TSMC improved w.r.t. the TSMC 82:4% and SDDRE 47:2%. Based on the presented data in Shiri et al., 12 for t f = 4 s, the proposed TSMC response has been validated. Setting the final time as an input with two sliding surfaces was proposed in Shiri et al., 12 though in this current work, a unified sliding surface is considered with additional sub-optimality in the gain.
The chattering effect on the input is always an issue in the implementation of different kinds of sliding mode controllers. Here we also observe the chattering in the input signals, Figure 4

FWFR
The characteristics of the robot define the behavior of the system; the weight of an FWFR prototype is usually less than 1 (kg). The aerodynamics force and drag in the FWFR case play an important role. Unlike aircraft, during the regulation in X axis, there are significant variations in other directions as well, and therefore, aerodynamics force and drag should be included in three axes. The tuning of the controller is sensitive and dependent on the initial and final conditions. The physical characteristics of a sample FWFR are reported in Table 2.
The actuator limitations are applied by À(p=6) ł d (Á) ł (p=6) and À0:5248 ł F T ł 7:8725. The elements of the mixer matrix are also considered as: Then, the inputs could be found by u = M À1 x U, where the limits are applied by The main effects of disturbance caused by flapping are in z c ; hence, the disturbance vector is considered in the system (equation (35)) to complete the modeling, in which v = 30 (rad=s) is the flapping frequency and d F = 0:1. Based on the assumption of the modeling, the trajectory should be long enough in X axis direction. The initial condition is regarded as zero except the initial velocity of u(0) = 1 m=s and angle u(0) =À (p=12) (rad), which means throwing the FWFR to the air, x(0) = ½0, 0, 0, 0, À p=12, 0, 1, 0, 0, 0, 0, 0 T . The final condition is x(t f ) = ½15, À 5, 2, 0:344, À 0:132, À 0:321, 0 1 3 6 T and the simulation time is chosen 15 (s).
The control gain selection is quite sensitive due to the flight conditions and limitations of the model. The tuning process includes SDRE and TSMC gain selection. The first step is to define the weighting matrices of the SDRE without considering the TSMC; then, with stable SDRE tuning, the TSMC must be   Tuning the TSMC also presents h = ½1, 1, 1, 0:1, 0:1, 0:05 T , b = 1, q = 7, and p = 9 (q and p are the power of absolute terms in equation (8)). For the combination of the TSMC and SDDRE, the orientation weighting matrices must be changed to Q o = I 6 3 6 . Simulating the model, the position variables are presented in Figure 7. The orientation variables of the system are depicted in Figure 8. The desired orientation angles were defined based on the geometry of the initial and final position to keep the system aligned in a straight line that connects both points. The linear and angular velocities of the FWFR are presented in Figures 9 and 10. The control input force and moments are presented in Figure 11. The nature of the TSMC imposes chattering on the control signals that could be amended by the use of saturation function; however, this reduces the precision in finite-time control. The trajectory of the system is shown in Figure 12. The rudder angles of the left and right wings are shown in Figure  13(a) and (b) with respect. The variation of the lifting rudder on the rear edge of the wing is illustrated in Figure 13(c). The necessary generated thrust by flapping is shown in Figure 13(d). Finally, the angles of upper/down and left/right rudders on the tail are presented in Figure 13(e) and (f). The error was found 0:4395 mm. The chattering removal of the inputs ( Figure 13) using tanh function and e = 0:1 presents practical inputs for rudders and thrusts ( Figure 14); with the expense of more error in regulation 1:8781 mm.

Comparison and validation of FWFR results
To check the performance of the proposed controller applied to complex systems, the flight simulation has been validated with other controllers as well. Similar to section ''Illustrative example,'' SDDRE represents an optimal finite-time controller and the conventional TSMC represents a robust one. The PD + gravity has also been selected to check the system with a simple structured controller. The gains of the controllers are similar to section ''FWFR'' for TSMC and proposed    TSMC. The PD + gravity controller gains are K P, t = 0:5 3 I 3 3 3 and K D, t = 1:5 3 I 3 3 3 for translation control; and K P, o = 0:005 3 I 3 3 3 and K D, o = 0:01 3 I 3 3 3 for the orientation part. The gains of the SDDRE were also set similar to section ''FWFR.'' Conducting the simulation with the same conditions in section ''FWFR,'' the results are presented in Figure 15. The error of the regulation by the mentioned controllers is presented in Figure 15. The details of the errors are reported in Table 3. The proposed method gained the least error 0:43 mm using sign function. The chattering-free version of the controller also gained 1:87 mm error which is acceptable concerning 15 m regulation distance. Since the time of regulation was long enough for all the controllers, the difference between the proposed TSMC 0:43 mm and conventional TSMC 0:53 mm is negligible. However, decreasing the final time will show better results for the    proposed TSMC similar to section ''Illustrative example.''

Conclusion
A novel TSMC approach using the SDDRE is introduced. Finite-time control has been always a good option in control engineering. The conventional TSMC is a continuous finite-time controller; however, the desired final time is not commanded directly. The SDDRE presents a finite control gain that penalizes the states near the final time. Instead of using the standard SDDRE control law, the gain is extracted and partitioned to be used in the definition of the sliding surface.
The new nonlinear sliding surface presented a novel TSMC with the possibility of finite-time control as an input parameter. The Lyapunov stability criteria were used to deliver a stable adaptive gain to the TSMC. The proposed controller was validated by the simulation of a van der Pol oscillator. 12 The results were similar to the available data and also proposed TSMC scored 61% improvement in comparison with conventional TSMC and 61:3% with SDDRE; for the 6 s final time.
To check the performance of the controller, a complex model was also simulated. The FWFR possesses nonlinear uncertain dynamics which is a challenging case in the control field. The proposed TSMC controlled the model successfully and was compared with SDDRE and conventional TSMC to validate the results. The superiority of the proposed TSMC was also shown by 18:3% improvement concerning conventional TSMC, 93:5% over SDDRE and 96:8% over PD + gravity.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the European Research Council as part of GRIFFIN ERC advanced Grant 2017, Action 788247, and by the European Commission H2020 Program AERIAL-CORE, contract no. 871479. Partial funding received by the Plan Andaluz de Investigacia´n, Desarrollo e Innovacioen (PAIDI) 2020 through the Project HOMing Pigeon bOT (HOMPOT) under Grant PY20_00597

ORCID iD
Saeed Rafee Nekoo https://orcid.org/0000-0003-1396-5082  The bold value represents the best results, obtained by the proposed method, a combination of the TSMC and SDDRE using the sign function.