Closed-loop nonlinear optimal control design for flapping-wing flying robot (1.6 m wingspan) in indoor confined space: Prototyping, modeling, simulation, and experiment

Please cite this article as: S.R. Nekoo and A. Ollero, Closed-loop nonlinear optimal control design for flapping-wing flying robot (1.6 m wingspan) in indoor confined space: Prototyping, modeling, simulation, and experiment. ISA Transactions (2023), https://doi.org/10.1016/j.isatra.2023.08.001. ISA Transactions xxx (xxxx) xxx Contents lists available at ScienceDirect ISA Transactions journal homepage: www.elsevier.com/locate/isatrans Practice article Closed-loop nonlinear optimal control design for flapping-wing flying robot (1.6 m wingspan) in indoor confined space: Prototyping, modeling, simulation, and experiment✩ Saeed Rafee Nekooa,∗, Anibal Olleroa,b aThe GRVC Robotics Lab., Departamento de Ingeniería de Sistemas y Automática, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Seville, 41092, Spain bFADA-CATEC, Centro Avanzado de Tecnologías Aeroespaciales, Seville, 41300, Spain article info Article history: Received 20 April 2023 Received in revised form 1 August 2023 Accepted 2 August 2023 Available online xxxx Keywords: Flapping-wing robot Aerial robotics SDRE Nonlinear Optimal Closed-loop control Flight control abstract The flapping-wing technology has emerged recently in the application of unmanned aerial robotics for autonomous flight, control, inspection, monitoring, and manipulation. Despite the advances in applications and outdoor manual flights (open-loop control), closed-loop control is yet to be investigated. This work presents a nonlinear optimal closed-loop control design via the state-dependent Riccati equation (SDRE) for a flapping-wing flying robot (FWFR). Considering that the dynamic modeling of the flapping-wing robot is complex, a proper model for the implementation of nonlinear control methods is demanded. This work proposes an alternative approach to deliver an equivalent dynamic for the translation of the system and a simplified model for orientation, to find equivalent dynamics for the whole system. The objective is to see the effect of flapping (periodic oscillation) on behavior through a simple model in simulation. Then the SDRE controller is applied to the derived model and implemented in simulations and experiments. The robot bird is a 1.6 m wingspan flapping-wing system (six-degree-of-freedom robot) with four actuators, three in the tail, and one as the flapping input. The underactuated system has been controlled successfully in position and orientation. The control loop is closed by the motion capture system in the indoor test bed where the experiments of flight have been successfully done. ©2023 The Author(s). Published by Elsevier Ltd on behalf of ISA. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Flapping-wing flying robots have been increasingly used in research and autonomous robotics targeting the applications of monitoring and inspection [1,2], manipulation and sampling [3], an alternative way of flight with less noise in nature [4], etc. One application and challenge in FWFR is perching on a branch. Perching is needed to perform a post-perching activity such as sampling, manipulation, taking a picture, or other tasks. Perching on a branch was done using multi-rotor drones [5], fixed-wing unmanned aerial vehicles [6], flapping-wing robot [7], etc. The flight and perching by the ornithopter used a series of linear controllers to regulate the height of the robot close to the branch position [7]. The flight was also performed in an indoor testbed limited to a 15 (m) flight path. Linear controllers are preferable for ✩Funding: This work was supported by the European Project GRIFFIN ERC Advanced Grant 2017, Action 788247. ∗Corresponding author. E-mail addresses: [email protected] (S.R. Nekoo), [email protected] (A. Ollero). control designs if they could do the task satisfactorily since they are fast in the computation of control loops. To move towards a model-based nonlinear control, an efficient representative model is demanded. The objective of this work is to control a robot bird in an indoor test bed. The flapping-wing robots are suitable prototypes for outdoor flights since they need a large area to freely maneuver. Control design with repeatability is difficult in outdoor tests due to a lack of precision in the measurement of position feedback. The indoor test bed conversely has the possibility of measurement and repetition of the flights using Opti-Track system. The indoor flight zone also offers an area with the same conditions (without wind and disturbance). The flapping-wing robot possesses a six-degree-of-freedom (DoF) model with four actuators. The actuators interact in all DoFs of the system, i.e. the tail could change the pitch angle, consequently, it affects the height. The flapping itself can ascend or descend the robot. Flapping is the source of forward thrust as well if the robot has a close-to-zero pitch angle. It is possible to control all the six-DoF variables of the system in outdoor unlimited space, though in this study, in indoor experiments, it is not in the interest of the work https://doi.org/10.1016/j.isatra.2023.08.001 0019-0578/©2023 The Author(s). Published by Elsevier Ltd on behalf of ISA. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/). S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx to consider longitudinal control since the flight zone is limited to 15 (m). This research has been done within the framework of the GRIFFIN project, ‘‘General compliant aerial Robotic manipulation system Integrating Fixed and Flapping wings to INcrease range and safety’’, advanced grant ERC. The GRIFFIN project targets the use of aerial robotics, specifically flapping-wing technology for building prototypes and robots for the application of monitoring and inspection, manipulation, and lightweight object transportation. The early-stage robot birds performed outdoor manual flight [8], gliding model, and flight [9]. The closed-loop control of the FWFR was done with linear designs and led to the possibility of perching on a branch autonomously [7]. This work will continue the research in this line of investigation by developing a more sophisticated control based on the state-dependent Riccati equation. The contribution of this work with respect to the previous research is the implementation of nonlinear optimal control in closed-loop design and the presentation of an equivalent dynamic to imitate the flapping effect on the six-DoF model of the robot. The modeling of flapping-wing robots is a harder challenge than the fixed-wing unmanned aerial platform due to the additional complexity through the flapping of wings which adds two more 3D objects to the main body of the robot and consequently the generated aerodynamics terms. This additional motion of the wings changes the dynamics from single-body to multibody system dynamics. Tu et al. investigated the flight control of flappingwing robots with damaged wings [10]. The presented dynamics included the drag and damping coefficients to consider the effect of wings and their motion. Guzman et al. studied the dynamics and aerodynamic coefficients of three types of tails for FWFRs; using two servomotor actuators with different configurations [11]. Fan et al. presented deformable wings for FWFRs and studied the effect of flexibility in the control performance [12]. A simplified model for the flapping wing was reported and validated experimentally, presenting the drag and lift terms due to actuation [13]. Despite the various versions of the dynamics of the robot birds, model-based controllers have been rarely implemented on the experimental platforms. The Lagrange method will be used in this work to derive an equivalent dynamic for translation and a conventional one for orientation. The equivalent dynamic will present the effect of periodic flapping on the base of the robot similar to experimental data. The derivation of the dynamics can be also done using other techniques such as Kane’s method [14,15], Newton–Euler approach [16,17], and Gibbs–Appell formulation [18–20]. Modeling the robotic systems using different dynamic methods is getting popular and embedding the recursive structure on them provides powerful and fast derivations for complex systems. The focus of this work is on the control of the FWFR. The beginning era of flapping wings presented many outdoor flights using radio controller (RC) transmitters. The term FWFR or ornithopter robots refers to unmanned systems (from a robotics point of view), which tried to perform the autonomous flight. The concept of flapping dates back to the ancient Greek legend of Daedalus and Icarus, Leonardo da Vinci, and the early years of flapping wing systems started in 1870 by Gustav Trouvé [21]. So it is worthwhile echoing the phrase of Chronister that the word ornithopter was defined already and several platforms accomplished the flight way back in history [21]. Flying outdoors even recently needed a human pilot due to the inaccuracy of the positioning systems and inertial measurement units (IMUs); Pan et al. presented two prototypes with a wingspan of more than 2 m flying outdoor [1,22]. Lee et al. designed and experimented with a robust controller for attitude control of a robot bird subjected to wind disturbance [23]. The outdoor control was also reported using global positioning system (GPS) feedback for localizing and an IMU for orientation [24]. The most popular reported controller for the flight was proportional–integral–derivative (PID) design [25], applied on trajectory tracking of a circular path [26], flappingwing miniature aerial vehicle [27], and other platforms; however, active disturbance rejection control [28], and adaptive control were also implemented on gliding mode of the FWFR [29]. Here in this work, the state-dependent Riccati equation controller is implemented on the attitude control of the flappingwing flying robot. The application of the SDRE controller diverged in the early 90s in aerospace control [30,31], then exploring other fields such as robotics [32,33], control of underwater platforms [34,35], etc. The SDRE and terminal sliding mode controller were also simulated on a version of a flapping wing robot though the experimental implementation was not exercised [36]. The SDRE is an optimal nonlinear closed-loop approach that needs a model for designing the controller. The available models for the flapping-wing systems are too complex for implementation and delivering a solution to the Riccati equation as an online onboard implementation on the robot bird. Here the SDRE is implemented on an equivalent dynamic whose output is similar to the experimental data of the FWFR in many (regulation) point- to-point control records. The experimental flight data of the robot will be also compared with the simulation results. Experimental implementation of nonlinear optimal controllers, specifically the state-dependent Riccati equation has been a challenge in practical works. From the early works in the SDRE experimental tests, one could refer to hardware-in-the-loop with two computers back in the late 90s [37], or control of an underactuated robot in 2002 [32]. In 2003, practical implementation of the SDRE was reported for unmanned helicopter control [38]. It was reported the Riccati was solved at each time-step onboard the flying system for tracking the commanded signal of the pilot (through RC transmitter) after takeoff. Here in this current work, pilot or RC transmitter is not used, and the control design is done on the onboard computer. The inverted pendulum was also an interesting experimental platform to show the capability of the SDRE controller [39]. The main challenge in the experimentation of the SDRE has been the solution to the Riccati at each timestep that led to several implementations of stationary setups [40]. Communication between the main processor and the robot onboard computer was an effective method to reduce the burden of the computational load from the experimental platform [41]. The time-step of simulation then increases in comparison with linear controllers; so far the reports showed an average of 0.05 to 0.085 (s) for rotary inverted pendulum [42], and variablepitch pendulum [43,44]. A PID-SDRE control design was studied for orientation control of a quadrotor drone [45]. The summary of this short discussion on the experimental implementation of the SDRE indicates that the beginning of practice started with stationary setups, communication between the computer and the processor of the platform, and then the implementation and solution to Riccati was done at each time-step onboard. Now the trend is to implement the algorithm totally onboard the platform to remove the risk of wireless communication, especially on flying systems. The contributions of this work are as follows: 1. The first major contribution is to implement the SDRE controller (theoretical and experimental study with onboard implementation) on the attitude control of the flappingwing flying robot for the first time using an equivalent dynamic model. The linear control on E-Flap was implemented in [7], and here a nonlinear closed-loop controller upgraded the system. To the best knowledge of the authors, the SDRE has been implemented experimentally on FWFR for the first time, in this work. 2 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx Fig. 1. (a): The CAD design of the flapping-wing flying robot in 3D view. (b): The zoomed view of the hinges, (c): and the tail motors. The elevators right and left are installed with an acute angle to avoid contact between the elevators in the actuation. 2. The second major contribution is the addition of Y-axis control to the height and orientation control, reported in previous works [7,25]. With this control design, the five- DoF of the system is controlled except for the X-axis which implies the necessity of forward flight to preserve the stability of the robot in regulation. Controlling the yaw angle of the robot results in an indirect control of the Y-axis; however, having control over that as an independent control variable would increase the possibility of more complex maneuvers. 3. The first minor contribution is the change in the design of the tail and adding one more actuator to have three servomotors for obtaining more control over the robot. It has been common to consider two servomotors for tail control [2,11,25,46,47]; and here for the first time a controlled flight is done with this specific tail design. The addition of actuators increases the weight and interaction of the control variables. A successful point-to-point control has been performed thanks to the new tail design using four actuators, flapping, rudder, and left/right elevators. 4. The second minor contribution is a technological one, using completely printed tails instead of conventional custommade tails with fabrics and carbon fiber (CF) rods as ribs. This will facilitate a lot the fabrication process, symmetry in the bird, and repair of the system in case of a crash. The paper structure: Section 2presents the mechanical design procedure and mechatronics of the flapping-wing robot. Section 3 describes the modeling of the robot bird and the idea of the equivalent dynamics in the flapping as a periodic excitation of the base. The controller design and the structure of the statedependent Riccati equation are reported in Section 4. Sections 5 and 6present simulation, experimental flight data, and results respectively. Concluding remarks are reported in Section 7. 2. Mechanical and mechatronics design 2.1. Mechanical design The objective of modification and development of the new prototypes is to deliver a more complete system and reach better performance while removing the drawbacks in previous versions — an iterative motion in the evolution of robot birds. This FWFR is built in the context of the GRIFFIN project,1as a new bird with the following changes with respect to E-Flap [25]: (1) increasing the wingspan from 1.5 m to 1.6 m to add lift force, (2) addition of one actuator to the tail to enhance maneuverability, (3) present a combination of ‘‘V’’ tail and ‘‘T’’ tail, please see Ref. [11], to achieve more control over the robot, (4) using the idea of printed tails instead of building the tails with fabrics and carbon fiber rods, (5) obtaining the Y-axis control and implementation of a nonlinear optimal controller, SDRE, onboard the robot. The 3D printed tails are more accurate, it is easy to replace them and it increases the symmetry of the design. Making the tails by fabric requires cutting the fabric based on a pattern and use of glue which puts some inaccuracy in manufacturing. Both ‘‘V’’ and ‘‘T’’ tails in previous prototypes were built with two servo motors while here three servomotors are considered in the design. The 3D CAD design of the FWFR is shown in Fig. 1. The design is done to consider two processes, CNC cutting of CF plates and 3D 1https://griffin-erc-advanced-grant.eu/ 3 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx Fig. 2. The side view of the CAD design of the flapping-wing flying robot. Fig. 3. (a): The wing structure of the robot, main tube, and the ribs, 1.6 (m) wingspan; (b): the old version using 3D printed part and glue; (c): the new version passing the ribs through the tube without using glue. Fig. 4. The tail structure of the robot, definition of the motions and torque and thrust reactions; (a) shows both elevators up which results in negative pitch angle and consequently upward thrust force; (b) placed both elevators down that generates positive pitch angle and negative thrust force; (c) sets left elevator up and the right one down to produce negative roll angle; (d) sets the left elevator down and the right one up to generate positive roll angle; (e) places the rudder to right that provides negative yaw angle; (f) sets the rudder to left that provides positive yaw. printing of parts with polylactide (PLA) material. The main plates for holding the motors and gears were made from CF plates of 1.5 (mm) thickness, presented in Fig. 2. The brushless DC motor is fixed on the left-side plate, gears were installed using bearings on both sides, and the battery is set between the two plates at the tip of the robot to balance the weight distribution and keep the center of gravity near the flapping point. The gear reduction ratio is n=Z2 Z1×Z4 Z3=42 where Z1=12,Z2=72,Z3= 8,Z4=56 and it reduces the high-speed brushless DC motor to a proper zone for flapping frequency. The reduced rotary motion of the gears will be transformed into a reciprocating motion by a cam mechanism. The aluminum hinges hold the main wing tubes, with 6 (mm) and 4 (mm) outer/inner diameter made of CF. The ribs are 1.5 (mm) CF rods passing through the holes on the main 6 (mm) CF tubes, Fig. 3. In E-Flap version [25], the ribs passed through a 3D printed holder and the holder was installed (glued) on the main wing tube. This new design consideration reduces the 3D-printed holders and reduces the weight of the wings. Clearly, it has more accuracy and symmetry since gluing is removed from this part of the manufacturing process. The position of the robot is measured using the Opti-Track system. Seven emitting balls are installed on the robot bird to reflect infrared light to the Opti-Track cameras with high frequency. Two markers will be installed at the tip of the FWFR on the blue part antenna shape, see Fig. 1. Another pair will be installed above the hinges, one above the processor board, and the last two beside the servomotors. The tail of the robot includes three servomotors for two elevators and one rudder. The combination of the three actuators contributes to roll, pitch, and yaw (orientation) control. The elevators are also effective in thrust production. The motion is briefly explained as follows. Placement of both elevators up results in a negative pitch angle and consequently upward thrust force, Fig. 4-(a). Then this actuation increases the height of the system in flight (increasing positive Z-axis motion). The opposite 4 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx motion regulates the robot bird down and reduces the height in flight by placing both elevators down, Fig. 4-(b). Setting the left elevator up and the right one down will generate a negative roll angle and rotates the bird around X-axis, Fig. 4-(c). Setting the left elevator down and the right one up generates a positive roll angle, Fig. 4-(d). The rudder has a role in yaw rotation. Placing the rudder to the right provides a negative yaw angle and setting the rudder to the left provides positive yaw, Figs. 4-(e) and -(f). It will be shown that the rudder will contribute to the Y-axis control in the modeling section. The last mechanical part is the linkage between the body and the tail which is constructed by a CF plate of 1.5 (mm) with 15 (mm) separation. In the whole design of the plates, extra material was removed by cutting different shapes and holes in them to reduce the mass as much as possible. 2.2. Mechatronics and electronics design The motion of the flapping wing is generated by a brushless DC motor, Hacker A20-26M kv1130. It has a direct feed from the driver ‘‘Tmotor F 35 A 32 bit 3S’’ which can work with 3S to 5S lipo batteries. The schematic view of the electronics and actuators of the system are presented in Fig. 5. Here the power supply is chosen a 4S1P with 450 (mAh) and 14.8 (V) operating voltage. The input control signal to the DC motor is received from a driver board, PCA9685, which is capable of handling up to 16 pulse-width modulation (PWM) outputs. The operating voltage of this module is 5V and it communicates with the Raspberry PI 4B with serial peripheral interface (SPI) protocol. The power supply to the PCA9685 is provided by LM7805 conventional regulator, with 5V voltage and up to 1.5(A) current. This output current will be consumed by three servomotors for the tail of the FWFR, SH-0255MG type. These small servomotors can work with operative voltage [4.8,6.0](V) with corresponding output torque [3.1,3.9](kg.cm) for the input voltage range. The actuation time in the range of [4.8,6](V) is [0.13,0.16](s), with operating travel π/2 (rad) from 1000 →2000 (µs). These servo motors do not provide position feedback for the control, and thanks to the external measurement system, Opti-Track, angular measurement of the tails is not necessary. The position and orientation of the robot bird will be obtained by the Opti-Track system with an accuracy of 1 to 2 (mm), depending on the calibration of the system. The processor board of the flapping-wing robot is a Raspberry PI 4B which is a conventional small single-board computer. It will be powered up by a voltage regulator XL4015 which can work with input voltage [8,36](V) and deliver the output voltage of [1.25,32](V) with a current of up to 5(A). It is important to set the input power of the Raspberry PI at 5.2 (V) to avoid the problem of restarting and low-voltage warning. The processor communicates wireless with the motion capture system through wireless fidelity (WiFi) for receiving the position and orientation feedback of the robot in flight for employing in the closed-loop control loop. Therefore, independent feedback from the actuators is not necessary for the controller which is crucial since brushless DC and servomotors do not usually provide position feedback. The motion capture system was installed in the testbed, confined space of 20 ×15 ×7 (m), with 28 cameras that provided approximately 15 (m) diagonal flight path for the ornithopter. The communication with the Opti-Track system is done by the robot operating system (ROS), or to be more accurate, ROS1, to read the position and orientation of the FWFR body. Virtual-reality peripheral network (VRPN) interface is used for reading the data of the motion capture system. It will be read by a Python script code in the control loop by subscribing to the internet protocol address (IP) of the robot bird defined by Motive software. 2.3. Mass distribution and details This section presents the mass breakdown of the components of the FWFR. The previous version of the robot bird weighed 500 (g), E-Flap prototype reported in [25]. The FWFR in this work is 567.3 (g) without the battery, shown in Fig. 6. Here we highlight the changes that led to the additional weight of this version. The processor of the E-Flap was Nano PI NEO which was 9.7 (g) changed to Raspberry PI 4B with the weight of 46 (g). The voltage regulators were installed on the main CF plate of the robot in E- Flap though here we used LM7805 and XL4015. The wingspan of the robot increased 10 (cm) and another servo has been added to the robot on the tail. These modifications were done to simplify the repair and manufacturing of the robot and add the payload by increasing the wingspan. The weight breakdown of the FWFR is presented in Fig. 7. The breakdown shows that the majority of the mass is devoted to the structure and mechanical parts, 42.5%, and in the second place, the wing with 20.9% share of the total weight. Using PWM module PCA9685 allows us to add servo inputs up to 16, which will be used for adding more servomotors for manipulation and leg of the robot birds in future works. 3. Modeling and equivalent dynamics The dynamics modeling of the flapping-wing flying robot has the additional complexity of flapping with respect to airplane dynamics. A complete six-DoF modeling is considered for the translation and orientation of the robot bird. The dynamics of the airplanes have been well-developed in the literature; however, the current valuable dynamics of flapping-wing robots are too complex to be used in model-based control designs. One approach is numerical model development for flapping wing dynamics and computational fluid dynamics (CFD) [48]. The complexity of the flapping-wing model is rooted in the definition and computation of aerodynamics terms [49]. To present a model for the implementation of a nonlinear controller, a simplified approach is chosen. In order to consider this approach, it should be noted that a flapping-wing robot is different from an airplane dynamics: A1. The ratio of the weight to the area of the wing is much lower than in an aircraft; in other words, an aircraft is heavier with respect to the surface of its wings. A2. The speed of the forward flight of an airplane is way higher than FWFR. A3. The FWFR has a periodic oscillation in the flight path synchronous with the frequency of the flapping. A4. The flight condition of this work is an indoor confined space, without wind and external disturbance. Moreover, the following assumptions are considered to facilitate the modeling process: B1. The robot flies in the forward direction, X(t), with forward velocity more than the minimum required speed to keep the stability of flight, ˙ X(t)≥˙ Xmin [36]. B2. There is always an initial speed for launching the robot which means an initial velocity ˙ X(0) ≥˙ Xmin [36]. B3. The aerodynamics effect is neglected in the modeling due to low-speed flight conditions. This will reduce the precision of the modeling though it is intentionally done to simplify the control design. B4. The effect of the flapping is considered by equivalent dynamics and base excitation modeling. The actuators of the FWFR are flapping wings and the motion of the tail (rudder and elevators). The main flapping provides the source of lift force to compensate for the gravity and the forward 5 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx Fig. 5. The electronics and actuators block diagram of the FWFR. Fig. 6. The weight (mass) measurement of the FWFR. Fig. 7. The weight (mass) breakdown of the FWFR; the unit of measurement is (g) in the chart. Fig. 8. The model of the flapping-wing robot, axes definition, and actuation. thrust force. The two elevators are defined by {δER(t), δEL(t)}(rad) for right and left one respectively. The motion of the rudder is denoted by δR(t) (rad), presented in Fig. 8. The combination of the motion of the tail components generates three main torques along the local axis (placed at the center-of-mass (CoM) of the robot), presented in Fig. 4. The tail also contributes indirectly to the change in height. By changing the pitch angle, the flapping will move the robot higher/lower during the forward motion. It is also true that the rudder can indirectly change the position of the robot in the Y-axis. The axes definition and the schematic of the robot are presented in Fig. 8. The global fixed coordinate is presented by {X,Y,Z}and it measures the local (moving) coordinate by ξ1(t)= [xc(t),yc(t),zc(t)]⊤(m). The orientation of the local axis is measured by ξ2(t)= [φ(t), θ(t), ψ(t)]⊤(rad) with respect to main fixed frame. The linear velocity of the moving axis is υ1(t)= [u(t), v(t), w(t)]⊤(m/s) and the angular velocity of that is υ2(t)= [p(t),q(t),r(t)]⊤(rad/s). The kinematic relationship between the global and local coordinates is defined by [50]: ˙ ξ1(t)=Rb(ξ2(t))υ1(t),(1) ˙ ξ2(t)=T(ξ2(t))υ2(t),(2) where Rb(ξ2(t)) =[cψcθ−cφsψ+sφsθcψsφsψ+cφsθcψ sψcθcφcψ+sφsθsψ−sφcψ+cφsθsψ −sθcθsφcθcφ], T(ξ2(t)) =[1sφtθcφtθ 0cφ−sφ 0sφ/cθcφ/cθ]. Neglecting the effect of aerodynamics, based on Assumption B3, the robot is considered a 3D flying object. The generalized coordinate vector for a flying object (the common form in the aerospace community) is selected as q(t)={xc(t),yc(t),zc(t), φ(t), θ(t), ψ(t)}. Applying the Lagrange method, the equations of motions are derived as: mb¨ ξ1(t)+mbge3=Rb(ξ2(t))FB(t),(3) J(ξ2(t))¨ ξ2(t)=τB(t)−c(ξ2(t),˙ ξ2(t)),(4) where mb(kg) is the mass of the base (total mass without the mass of wings; the weight of the wings will be added later in the base excitation model), g=9.81 (m/s2) is the gravity 6 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx Fig. 9. A sample of flight data for height regulation of the FWFR, experimental data for three set-points [7]; bold lines show the average motion. constant, τB(t)= [τφ(t), τθ(t), τψ(t)]⊤(N.m) is the counteracting applied torque vector to the rotational dynamics, FB(t)= [Fx(t),Fy(t),Fz(t)]⊤(N) is the counteracting applied force vector to the translation dynamics, both set on the local axis and e3= [0,0,1]⊤.J(ξ2(t)) :R3→R3×3is the inertia matrix of the rotational dynamics and c(ξ2(t),˙ ξ2(t)) :R3×R3→R3includes the Coriolis and centrifugal terms. Based on Assumption B1, the robot always flies at a forward speed in the X-axis and here the input control component (which would control the Xposition) of that is equal to zero, Fx(t)=0. In reality, the frequency of the flapping in a long-term flight will define the forward speed of the X-axis motion; however, in this 15 (m) short distance (the diagonal length of indoor confined testbed), the forward speed was always observed from experimental data u(t)= [2,6](m/s). Since the system is underactuated, the frequency of flapping was devoted to height control. For the sake of simplicity, a constant average speed of forward flight is considered for the modeling and simulation section. This average speed is based on the observation of experimental data from many recorded flights. Observing the experimental data, especially the log record of the Z-axis, presented in Fig. 9, it can be concluded that the height component of the generalized coordinate zc(t) is subjected to a periodic and persistent oscillation, with the frequency of [3.5,4.5](Hz) and amplitude of [10,15](cm). It was already mentioned that the modeling of flapping action in the dynamic equation of the FWFR is difficult, and the oscillation in height control brings up the idea of modeling the system with equivalent dynamic and base excitation [51], Fig. 10. The motion of the wing plays the role of oscillatory part of the dynamics and the body of the robot bird receives the forced excitation. Assumption 1. The frequency of the flapping is limited to [3.5,4.5](Hz) which is sufficient to change the altitude of the flapping wing robot [52]. The pitch angle of the FWFR varies between [0.2,0.5](rad) which is negligible though it is used to regulate the pitch angle slightly around an equilibrium point which provides flight stability for the bird. With a proper pitch angle around 0.4 (rad), the frequency of flapping ascends or descends the FWFR. The interactions of the motion of the elevators and flapping to height control are undeniable though devoting an independent role of height control to flapping simplifies the design. Assumption 2. The robot flies with negligible change in the pitch angle, then the role of flapping frequency is dominant in the change of the altitude [52]. Fig. 10. The equivalent dynamic of flapping action as a base-excitation system for the FWFR. Considering Assumptions 1 and 2, the system can be modeled as an equivalent dynamic with base excitation in which the oscillation of the flapping directly disturbs z-axis and indirectly, other degrees of freedoms through the nonlinear coupling of dynamic. The free-body diagram of the simplified model is depicted in Fig. 10. Regarding the balance between the forces in the z-axis, the equation of motion can be derived, more details could be revisited in [53,54]. Equation of motion of free-body diagram of Fig. 10, as an independent one-DoF mass–spring–damper system, is presented by: mb¨ zc(t)+c˙ zc(t)+kzc(t)+(mb+mw)g =mw¨ zw(t)+c˙ zw(t)+kzw(t)+Fz(t),(5) where mwis the mass of wings, c(N s/m) is the damping coefficient, k(N/m) is the elastic constant, zw(t)∈R,(m) is the base excitation motion. mb¨ zc(t) and mw¨ zw(t) are the effect of the acceleration on the mass of base and wings. In the freebody-diagram, c(˙ zc(t)−˙ zw(t)) is the damping force and k(zc(t)− zw(t)) the generated force due to elastic characteristics of the connection between the wings and the body. (mb+mw)gis caused by the gravity and Fz(t) is the lift force that compensates for the weight. The idea of equivalent dynamics and base excitation is expressed for a one-DoF system, only z-axis; however, this periodic disturbance will be incorporated into six-DoF dynamics which later will show the effect of excitation on other degrees of freedom. The base excitation is a periodic oscillation that can be presented by: zw(t)=z0sin (ω(t)t), ˙ zw(t)=z0{ω(t)+ ˙ω(t)t}cos (ω(t)t), ¨ zw(t)=z0[{2˙ω(t)+ ¨ω(t)t}cos (ω(t)t)−{ω(t)+ ˙ω(t)t}2sin (ω(t)t)], (6) in which z0(m) is the amplitude of the base excitation, and ω(t)∈ R,(rad/s) is the flapping frequency. If one divides (5) by mband considers damping ratio and natural frequency as ξ=c 2mbωnand ωn=√k mb(rad/s), respectively, it can be rewritten as: ¨ zc(t)+2ξωn˙ zc(t)+ω2 nzc(t)= z0([α{2˙ω(t)+ ¨ω(t)t}+2ξωn{ω(t)+ ˙ω(t)t}]cos (ω(t)t)+ [ω2 n−α{ω(t)+ ˙ω(t)t}2]sin (ω(t)t)) +Fz(t) mb−(1 +α)g,(7) which presents the standard form of second-order mass–spring– damper dynamic with forced input. Considering the effect of (5) 7 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx and (6) in the dynamic of six-DoF model, one can rewrite (3) in the form of: mb¨ ξ1(t)+mbge3=Rb(ξ2(t))[FB(t)+Ew(zc(t),zw(t),t)],(8) where Ew(zc(t),zw(t),t)=⎡ ⎣ 0 0 mw¨ zw(t)+c˙ zw(t)+kzw(t)−c˙ zc(t)−kzc(t)⎤ ⎦. It should be noted that gravity and external force are already present in (8) and the effect of base excitation will apply only in the Z-axis direction. It is true that due to the interaction of the rotation matrix in the translation dynamics, a smaller scale of oscillation will appear in the other X- and Y-axis as well. Selecting the state vector of the system as x(t)= [ξ⊤ 1(t),ξ⊤ 2(t),υ⊤ 1(t),υ⊤ 2(t)]⊤,(9) and derivative of x(t), result in the state-space representation of the system: ˙ x(t)=⎡ ⎢ ⎢ ⎢ ⎣ ˙ ξ1(t) ˙ ξ2(t) ˙ υ1(t) ˙ υ2(t) ⎤ ⎥ ⎥ ⎥ ⎦ .(10) Computing the time derivative of (1) and (2), and extracting ˙ υ1(t) and ˙ υ2(t) provide: ˙ υ1(t)=R−1 b(ξ2(t)) (¨ ξ1(t)−˙ Rb(ξ2(t))υ1(t)),(11) ˙ υ2(t)=T−1(ξ2(t)) (¨ ξ2(t)−˙ T(ξ2(t))υ2(t)).(12) Substituting (3),(4),(11), and (12) into (10) represents the state-space of the FWFR as (the argument tis not written for the brevity of the representation; the effect of base excitation is also considered in the state-space form): ˙ x=f(x)+g(x,u) ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Rb(ξ2)υ1 T(ξ2)υ2 R−1 b(ξ2)(1 mbRb(ξ2)Ew(zc,zw)−mbge3−˙ Rb(ξ2)υ1) T−1(ξ2)(−J−1(ξ2)c(ξ2,˙ ξ2)−˙ T(ξ2)υ2) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦+ ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 1 mbFB T−1(ξ2)J−1(ξ2)τB ⎤ ⎥ ⎥ ⎥ ⎦ . (13) 4. Controller design: The state-dependent Riccati equation 4.1. Nonlinear sub-optimal control design for orientation dynamics Consider a time-invariant nonlinear affine-in-control system: ˙ x(t)=A(x(t))x(t)+B(x(t))u(t),(14) where the input vector of the system is u(t)∈Rm, the state variables are collected in vector x(t)∈Rn, the state-dependent coefficient (SDC) parameterization of the system is represented by A(x(t)) :Rn→Rn×nand B(x(t)) :Rn→Rn×m.g(x(t),u(t)) = B(x(t))u(t) and f(x(t)) =A(x(t))x(t)in(14) are vector-valued smooth functions, piecewise continuous, that satisfies a local Lipschitz condition and: f(x(t)) :Rn→Rn, g(x(t),u(t)) :Rn×Rm→Rn. The objective of optimal control is to minimize the cost functional integral in an infinite-time horizon: J(·)=1 2∫∞ 0{u⊤(t)R(x(t))u(t)+x⊤(t)Q(x(t))x(t)}dt,(15) where Q(x(t)) :Rn→Rn×nand R(x(t)) :Rn→Rm×mare positive-semidefinite and positive-definite symmetric weighting matrices, respectively. More details on controllability and observability could be revisited in [55]. Constructing the Hamiltonian function: H(x(t),u(t),λ(t)) =λ⊤(t){A(x(t))x(t)+B(x(t))u(t)}+ 1 2[u⊤(t)R(x(t))u(t)+x⊤(t)Q(x(t))x(t)], and stationary condition of optimal control, ∂H(x(t),u(t),λ(t)) ∂u(t)=0, result in u(t)= −R−1(x(t))B⊤(x(t))λ(t).(16) Defining the co-state vector as λ(t)=P(x(t))x(t) and substituting that in (16), the SDRE control law is found [55]: u(t)= −R−1(x(t))B⊤(x(t))P(x(t))x(t),(17) where P(x(t)) :Rn→Rn×nis the solution to the SDRE and K(x(t)) =R−1(x(t))B⊤(x(t))P(x(t)) is the sub-optimal gain of the controller. The optimality condition ∂H(x(t),u(t),λ(t)) ∂λ(t)=˙ x(t) confirms (14), and the optimality condition ∂H(x(t),u(t),λ(t)) ∂x(t)= −˙ λ(t) provides: 1 2u⊤∂R(x) ∂xu+Q(x)x+1 2x⊤∂Q(x) ∂xx+A⊤(x)λ+ x⊤(∂A(x) ∂x)⊤λ+u⊤(∂B(x) ∂x)⊤λ= −˙ P(x)x−P(x)˙ x. (18) Substituting control law (17), system (14) and co-state vector λ(t) into Eq. (18) and mathematical manipulation, result in: 1 2(∂R(x) ∂xR−1(x)B⊤(x)P(x)x)⊤R−1(x)B⊤(x)P(x)x+Q(x)x+ 1 2(∂Q(x) ∂xx)⊤x+A⊤(x)P(x)x+(∂A(x) ∂xx)⊤P(x)x+ (∂B(x) ∂xR−1(x)B⊤(x)P(x)x)⊤P(x)x+ ˙ P(x)x−P(x)B(x)R−1(x)B⊤(x)P(x)x+P(x)A(x)x=0. (19) Eq. (19) provides state-dependent Riccati equation: A⊤(x)P(x)−P(x)B(x)R−1(x)B⊤(x)P(x)+Q(x)+P(x)A(x)=0, (20) and ‘‘the necessary condition for optimality’’ [56]: 1 2(∂R(x) ∂xR−1(x)B⊤(x)P(x)x)⊤R−1(x)B⊤(x)P(x)+ 1 2(∂Q(x) ∂xx)⊤+(∂A(x) ∂xx)⊤P(x)+ (∂B(x) ∂xR−1(x)B⊤(x)P(x)x)⊤P(x)+˙ P(x)=0, which collects the derivative terms of the result of Hamiltonian derivation ∂H(x(t),u(t),λ(t)) ∂x(t)= −˙ λ(t). Finite-time or infinite-time horizon control: The proposed controller for practical implementation is chosen infinite horizon SDRE. The infinite-time design and algebraic Riccati equation are based on tf→ ∞ in the cost function integral (15); and results 8 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx in, ˙ P(x)=0. It is true that for a short limited time, finite-horizon SDRE could be used, though it needs a more complex solution to the differential Riccati equation. In this work, for the sake of obtaining a shorter sampling time, an infinite-time horizon version is used. The SDRE (20) generates a solution to Riccati, a positivedefinite sub-optimal matrix for the gain of the control law (17). The SDRE controller is used for the orientation control of the FWFR. So, if the orientation dynamic of the FWFR is collected into a subsystem, presented as: ˙ x=[T(ξ2)υ2 T−1(ξ2)(−J−1(ξ2)c(ξ2,˙ ξ2)−˙ T(ξ2)υ2)]+[0 T−1(ξ2)J−1(ξ2)τB], hence, the SDC parameterization of the orientation dynamics is designed as: A(x)=[03×3T(ξ2) 03×3−T−1(ξ2)J−1(ξ2)C(ξ2,˙ ξ2)T(ξ2)−T−1(ξ2)˙ T(ξ2)], (21) B(x)=[03×3 T−1(ξ2)J−1(ξ2)],(22) where c(ξ2(t),˙ ξ2(t)) =C(ξ2(t),˙ ξ2(t))˙ ξ2(t) in which C(ξ2(t),˙ ξ2(t)) :R3×R3→R3×3. This process is also referred to as apparent linearization which is a factorization of the state vector x(t) from f(x(t)) and the input vector u(t) from g(x(t),u(t)). Remark 1. The Coriolis term −T−1(ξ2)J−1(ξ2)C(ξ2,˙ ξ2)T(ξ2) in SDC matrix (21), is multiplied by T(ξ2) from right-hand side which is not presented in state-space form (13). If ˙ ξ2is factored from Coriolis vector c(ξ2,˙ ξ2), it can be seen that C(ξ2,˙ ξ2)˙ ξ2is not compatible with state-vector (9). So, the kinematic Eq. (2) is substituted in that to form C(ξ2,˙ ξ2)T(ξ2)υ2. Remark 2. The derivative of rotation matrices, ˙ T(ξ2(t)) and ˙ Rb(ξ2(t)) are numerically computed in the control loop of the simulations. The control law of the orientation dynamic is presented as τB(t)= −R−1(x(t))B⊤(x(t))P(x(t))e(t),(23) where e(t)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ φ(t)−φdes θ(t)−θdes ψ(t)−ψdes ˙ φ(t)−˙ φdes ˙ θ(t)−˙ θdes ˙ ψ(t)−˙ ψdes ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,(24) shifts the equilibrium point from zero to a constant value in the workspace of the system; subscript ‘‘des’’ stands for the desired value. The gain P(x(t)) in (23) is a solution to the SDRE (20), updated by SDC (21),(22), weighting matrices R(x(t)) and Q(x(t)) of appropriate dimension. The block diagram of the control design for simulation is presented in Fig. 11. It demonstrates the orientation control using the SDRE which is fed to the orientation dynamic and the PD and PI control of the translation dynamic which will be explained in Section 4.2. 4.2. Linear control design for translation dynamic An average forward speed is considered for the X-axis motion in control design and simulation, based on the flight assumptions in Section 2. In the experiments, it happens in reality like that; in a short 15 (m) flight zone with an initial launching speed of 4 (m/s), the speed of the robot in X-axis reduces slightly from 4 (m/s) and stays around that value. Therefore Fx(t)=0 is set to zero. The error of the Z-axis is set as ez(t)=zc(t)−zdes and the one for Y-axis is defined by ey(t)=yc(t)−ydes. A proportionalderivative (PD) controller is considered for lateral control: Fy(t)= −KPey(t)−KD˙ ey(t),(25) where KP∈R+and KD∈R+are proportional and derivative gains of PD control. Considering Assumptions 1 and 2,ωmin ≤ω(t)≤ωmax in which ωmin =3.5(Hz) is the minimum frequency of the flapping, ωmax =4.5(Hz) is the maximum one, and ω(t) is the flapping frequency which is one of the control inputs. This limited change in the frequency of the oscillation ascends and descends the FWFR which means the elastic term and gravity are compensated indirectly by input lift force produced by flapping [52]: Fz(t):= mbgR−1 b(ξ2(t))e3+kzc(t)+τ(t),(26) where τ(t)=(ω(t)−ωmin)(pmax −pmin) (ωmax −ωmin)+pmin , that maps ω(t) from range of [21.9911,28.2743](rad/s) to the range of p∈ [−0.5,0.5]. A proportional–integral (PI) control finally regulates the frequency around a local equilibrium point for the flapping frequency ω0=4(Hz) =25.1327 (rad/s) which keeps the robot flying steadily in the air: ω(t)= −KP,fez(t)−KI,f∫t 0 ez(σ) dσ+ω0,(27) where KP,f∈R+and KI,f∈R+are proportional and integral gains of PI control. A PD controller has been chosen for the lateral control, the servomotor of the rudder. Since the design is done for a short final time of 4 (s), simple control structures are favorable to simplify the tuning. The tuning of the PD control was done through a series of iterations. The height control was sensitive to derivative terms and presented acceptable results with PI control. One of the difficulties in the control design was the fact that within the 15 (m) flight distance, the robot will be still inside the transient zone of the control design, which will be discussed more in Section 5. The robustness of the controllers, SDRE, PD, and PI was not studied in this work and the SDRE is in nature a nonlinear optimal method. So, the advantage of this approach is in optimality rather than robustness. The proposed control design is a combination of the SDRE for orientation and linear controllers for translation, Y- and Z-axis. Both controllers are popular methods, applied on a variety of systems, i.e. the application of SDRE could be found in [55]. This specific design showed a successful performance in the simulations which will be verified in the experiments as well. 5. Simulation results The time of the simulation is considered 4 (s) similar to the real results, the initial condition is set as x(0) =[−4.7,3.5,1.8, 0,0.15,−0.67,4,01×5]⊤and the final condition is chosen as x(4) = [9,−3.5,2,0,−0.55,0,01×6]⊤. The initial and final conditions of the flight are defined based on the experimental data and the confined space of the Opti-Track testbed to compare the model with the real data and prototype. The forward velocity of the robot in the X-axis is also set as u(t)=3.6 (m/s). The parameters of the robot are presented in Table 1. 9 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx Fig. 26. The snapshot of the flight experiment; the corresponding video of the flight is available as supplementary material of the paper. Fig. 27. The flight data record of 11 trials to show the repeatability of the control and implementation. More flights were performed and presented in the video presentation of the experiment, as supplementary material for the paper. law and system equation result in ˙ V(x(t)) = −x⊤(t)[Q(x(t)) + P(x(t))B(x(t))R−1(x(t))B⊤(x(t))P(x(t))]x(t) which is negative definite [59,60]. The same indoor condition provided the possibility of the repetitions of the experiments for more than 10 flights with statistical analysis, presented in Table 2. 7. Conclusions This work presented a theoretical and experimental study on the closed-loop control of a flapping-wing flying robot using a nonlinear optimal method, a state-dependent Riccati equation for orientation dynamic. The translation control was done through PD/PI control and the idea of equivalent dynamics for replicating the oscillation of the robot caused by flapping wings. The equivalent dynamic resulted in similar behavior of the robot with a simple and basic mathematical model to avoid complex machinery of aerodynamics. The mechanical design of the FWFR has been updated to increase the number of actuators, and an additional motor in the tail increased the controllability of the system. This change provided the possibility of the Y-axis control in addition to other degrees of freedom. The current controller design can regulate the system in 5-DoF due to this new change. The attitude control was done through a nonlinear control, which is a contribution with respect to previous works [25]. An analysis of the transient flight and equivalent dynamics for presenting a model for the FWFR was presented. The simulation and experimental results were compared and met a satisfactory alignment. The repetition of flights revealed an average error of 40 (cm) which is 2% of the flight trajectory. The video demonstration of the flights can be found as supplementary material of this work in the online version of the paper. Occasional unstable response due to a fault in tuning or the mechanical parts of the robot might lead to a crash; an example can be seen in the supplementary video that shows the deviation of the FWFR to the left while the rudder tries to correct the position to the right. The reason for that failure was two broken ribs of the left wing which caused a lack of lift of the left wing. Proposal for future works: The mechanical design of the robot could be refined to increase the mechanical resistance to possible impacts and to increase the performance of the flight. The mechanical enhancement is always limited by the weight of the components which imposes the system to be delicate. The flexibility of the wing and considering an elbow in a controlled manner also would help to enhance the lifting power and present a system closer to the real birds. From the control point-of-view, going towards six-DoF control in experimentation and real flight is a goal for future works. Declaration of competing interest The authors declare that there is no conflict of interest for this paper. Acknowledgments This work was supported by the European Project GRIFFIN ERC Advanced Grant 2017, Action 788247. The authors would like to thank Mario Hernandez for his help and useful discussions on the printing tail idea and manufacturing of the robot; Álvaro Satué Crespo for his help in electronics, Ernesto Sanchez Laulhe in aerodynamics; Raul Tapia and Juan Pablo Rodríguez Gómez for helping with the Opti-Track; and GRVC group members. The authors would like to thank the editors and reviewers of this paper as well. Appendix A. Supplementary data Supplementary material related to this article can be found online at https://doi.org/10.1016/j.isatra.2023.08.001. References [1] Pan E, Liang X, Xu W. Development of vision stabilizing system for a large-scale flapping-wing robotic bird. IEEE Sens J 2020;20(14):8017–28. [2] Wu X, He W, Wang Q, Meng T, He X, Fu Q. A long-endurance flapping-wing robot based on mass distribution and energy consume method. IEEE Trans Ind Electron 2022. 16 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx [3] Nekoo SR, Feliu-Talegon D, Tapia R, Satue AC, Martínez-de Dios JR, Ollero A. A 94.1 g scissors-type dual-arm cooperative manipulator for plant sampling by an ornithopter using a vision detection system. Robotica 2023;1–18. [4] Debiasi M, Lu Z, Nguyen QV, Chan WL. Low-noise flapping wings with tensed membrane. AIAA J 2020;58(6):2388–97. [5] Roderick WR, Cutkosky MR, Lentink D. Bird-inspired dynamic grasping and perching in arboreal environments. Science Robotics 2021;6(61):eabj7562. [6] Stewart W, Ajanic E, Müller M, Floreano D. How to swoop and grasp like a bird with a passive claw for a high-speed grasping. IEEE/ASME Trans Mechatronics 2022;27(5):3527–35. [7] Zufferey R, Tormo-Barbero J, Feliu-Talegón D, Nekoo SR, Acosta JÁ, Ollero A. How ornithopters can perch autonomously on a branch. Nature Commun 2022;13(1):7713. [8] Eguíluz AG, Rodríguez-Gómez J, Paneque J, Grau P, de Dios JM, Ollero A. Towards flapping wing robot visual perception: Opportunities and challenges. In: 2019 Workshop on research, education and development of unmanned aerial systems. IEEE; 2019, p. 335–43. [9] Lopez-Lopez R, Perez-Sanchez V, Ramon-Soria P, Martín-Alcántara A, Fernandez-Feria R, Arrue BC, et al. A linearized model for an ornithopter in gliding flight: Experiments and simulations. In: 2020 IEEE international conference on robotics and automation. IEEE; 2020, p. 7008–14. [10] Tu Z, Fei F, Liu L, Zhou Y, Deng X. Flying with damaged wings: The effect on flight capacity and bio-inspired coping strategies of a flapping wing robot. IEEE Robot Autom Lett 2021;6(2):2114–21. [11] Guzmán M, Páez CR, Maldonado FJ, Zufferey R, Tormo-Barbero J, Acosta JÁ, et al. Design and comparison of tails for bird-scale flapping-wing robots. In: 2021 IEEE/RSJ international conference on intelligent robots and systems. IEEE; 2021, p. 6358–65. [12] Fan X, Breuer K, Vejdani H. Wing fold and twist greatly improves flight efficiency for bat-scale flapping wing robots. In: 2021 IEEE/RSJ international conference on intelligent robots and systems. IEEE; 2021, p. 7391–7. [13] Sanchez-Laulhe E, Fernandez-Feria R, Ollero A. Simplified model for forward-flight transitions of a bio-inspired unmanned aerial vehicle. Aerospace 2022;9(10):617. [14] Banerjee A. Flexible multibody dynamics: algorithms based on Kane’s method. CRC Press; 2022. [15] Liu J, Chen X, Huang H, Ji S, Tu Q. A simplified method to analyze dynamic response of VLFS based on the kane method. J Offshore Mech Arct Eng 2023;145(3):031701. [16] Porez M, Boyer F, Belkhiri A. A hybrid dynamic model for bio-inspired soft robots—Application to a flapping-wing micro air vehicle. In: 2014 IEEE international conference on robotics and automation. ICRA, IEEE; 2014, p. 3556–63. [17] Bena RM, Yang X, Calderon AA, Perez-Arancibia NO. High-performance six- DOF flight control of the bee++: An inclined-stroke-plane approach. IEEE Trans Robot 2023;39(2):1668–84. [18] Zahedi A, Shafei A, Shamsi M. On the dynamics of multi-closed-chain robotic mechanisms. Int J Non-Linear Mech 2022;147:104241. [19] Honein TE, O’Reilly OM. On the gibbs–appell equations for the dynamics of rigid bodies. J Appl Mech 2021;88:074501. [20] Shafei A, Riahi M. The effects of mode shapes on the temporal response of flexible closed-loop linkages under the impulse excitation. Mech Syst Signal Process 2022;178:109256. [21] Chronister N. The ornithopter design manual, fifth edition. 2008, Published by the Ornithopter Zone, www.ornithopter.org. [22] Pan E, Xu H, Yuan H, Peng J, Xu W. HIT-hawk and HIT-phoenix: Two kinds of flapping-wing flying robotic birds with wingspans beyond 2 meters. Biomim Intell Robotics 2021;1:100002. [23] Lee J, Ryu S, Kim HJ. Stable flight of a flapping-wing micro air vehicle under wind disturbance. IEEE Robot Autom Lett 2020;5(4):5685–92. [24] Pan E, Liu J, Chen L, Xu W. The embedded on-board controller and ground monitoring system of a flapping-wing aerial vehicle. In: 2018 IEEE international conference on real-time computing and robotics. IEEE; 2018, p. 72–7. [25] Zufferey R, Tormo-Barbero J, Guzmán MM, Maldonado FJ, Sanchez- Laulhe E, Grau P, et al. Design of the high-payload flapping wing robot e-flap. IEEE Robot Autom Lett 2021;6(2):3097–104. [26] Wenfu X, Erzhen P, Juntao L, Yihong L, Han Y. Flight control of a largescale flapping-wing flying robotic bird: System development and flight experiment. Chin J Aeronaut 2022;35(2):235–49. [27] Roberts L, Bruck HA, Gupta SK. Autonomous loitering control for a flapping wing miniature aerial vehicle with independent wing control. In: International design engineering technical conferences and computers and information in engineering conference. 46360, American Society of Mechanical Engineers; 2014, V05AT08A013. [28] Liang S, Song B, Xuan J, Li Y. Active disturbance rejection attitude control for the dove flapping wing micro air vehicle in intermittent flapping and gliding flight. Int J Micro Air Veh 2020;12:1756829320943085. [29] Maldonado FJ, Acosta JÁ, Tormo-Barbero J, Grau P, Guzmán M, Ollero A. Adaptive nonlinear control for perching of a bioinspired ornithopter. In: 2020 IEEE/RSJ international conference on intelligent robots and systems. IEEE; 2020, p. 1385–90. [30] Mracek C, Cloutier J, Cloutier J, Mracek C. Full envelope missile longitudinal autopilot design using the state-dependent riccati equation method. In: Guidance, navigation, and control conference. 1997, p. 3767. [31] Mracek CP, Cloutier JR. Missile longitudinal autopilot design using the state-dependent riccati equation method. In: Proceedings of the international conference on nonlinear problems in aviation and aerospace. 1996, p. 387–96. [32] Erdem EB, Alleyne AG. Experimental real-time SDRE control of an underactuated robot. In: Proceedings of the 40th IEEE conference on decision and control (Cat. No. 01CH37228), vol. 3. IEEE; 2001, p. 2986–91. [33] Nekoo SR. Nonlinear closed loop optimal control: A modified statedependent riccati equation. ISA Trans 2013;52(2):285–90. [34] Naik MS, Singh SN. State-dependent riccati equation-based robust dive plane control of auv with control constraints. Ocean Eng 2007;34(11–12):1711–23. [35] Gao D-x, Lin X-h, Wang Y, Yang Q. Application of robust optimal sliding mode control for underactuated AUV system with uncertainties. In: 2020 Chinese control and decision conference. IEEE; 2020, p. 245–50. [36] Nekoo SR, Acosta J, Ollero A. Combination of terminal sliding mode and finite-time state-dependent riccati equation: Flapping-wing flying robot control. Proc Inst Mech Eng I 2023;237(5):870–87. [37] Langson W, Alleyne A. A stability result with application to nonlinear regulation: Theory and experiments. In: Proceedings of the 1999 American control conference (Cat. No. 99CH36251), vol. 5. IEEE; 1999, p. 3051–6. [38] Bogdanov A, Carlsson M, Harvey G, Hunt J, Kieburtz D, van der Merwe R, et al. State-dependent riccati equation control of a small unmanned helicopter. In: AIAA Guidance, navigation, and control conference and exhibit. 2003, p. 5672. [39] Suzuki S, Furuta K, Sugiki A, Hatakeyama S. Nonlinear optimal internal forces control and application to swing-up and stabilization of pendulum. J Dyn Sys Meas Control 2004;126(3):568–73. [40] Guzey UM, Copur EH, Ozcan S, Arican AC, Kocagil BM, Salamci MU. Experiment of sliding mode control with nonlinear sliding surface design for a 3-DOF helicopter model. In: 2019 XXVII International conference on information, communication and automation technologies. IEEE; 2019, p. 1–6. [41] Korayem MH, Lademakhi NY. Integrated nonlinear suboptimal control- and-estimator based on the state-dependent differential riccati equation approach. Optim Control Appl Methods 2022. [42] Nekoo SR. Digital implementation of a continuous-time nonlinear optimal controller: An experimental study with real-time computations. ISA Trans 2020;101:346–57. [43] Nekoo SR, Acosta JÁ, Heredia G, Ollero A. A benchmark mechatronics platform to assess the inspection around pipes with variable pitch quadrotor for industrial sites. Mechatronics 2021;79:102641. [44] Nekoo SR, Acosta JÁ, Heredia G, Ollero A. A PD-type state-dependent riccati equation with iterative learning augmentation for mechanical systems. IEEE/CAA J Autom Sin 2022;9(8):1499–511. [45] Giernacki W, Stepien S, Chodnicki M, Wróblewska A. Hybrid quasi-optimal PID-SDRE quadrotor control. Energies 2022;15(12):4312. [46] Jiang H, Zhou C, Xie P. Design and kinematic analysis of seagull inspired flapping wing robot. In: 2016 IEEE International Conference on Information and Automation. IEEE; 2016, p. 1382–6. [47] Ruiz C, Acosta JÁ, Ollero A. Optimal elastic wing for flapping-wing robots through passive morphing. IEEE Robot Autom Lett 2022;8(2):608–15. [48] Nakata T, Liu H, Bomphrey RJ. A CFD-informed quasi-steady model of flapping-wing aerodynamics. J. Fluid Mech 2015;783:323–43. [49] Han J-S, Breitsamter C. Aerodynamic investigation on shifted-back vertical stroke plane of flapping wing in forward flight. Bioinspiration Biomim 2021;16(6):064001. [50] Nekoo SR, Acosta J, Gomez-Tamm AE, Ollero A. Optimized thrust allocation of variable-pitch propellers quadrotor control: A comparative study on flip maneuver. In: 2019 Workshop on research, education and development of unmanned aerial systems. IEEE; 2019, p. 86–95. [51] Nekoo SR, Ollero A. Equivalent vertical dynamics of flapping-wing flying robot in regulation control: Displacement transmissibility ratio. In: 2023 International conference on unmanned aircraft systems. IEEE; 2023, p. 1301–7. [52] Nekoo SR, Ollero A. A proportional closed-loop control for equivalent vertical dynamics of flapping-wing flying robot. In: 2023 International conference on unmanned aircraft systems. IEEE; 2023, p. 1294–300. [53] Kelly JM, Konstantinidis D. Mechanics of rubber bearings for seismic and vibration isolation. John Wiley & Sons; 2011. [54] Sinha A. Vibration of mechanical systems. Cambridge University Press; 2010. [55] Nekoo SR. Tutorial and review on the state-dependent riccati equation. J Appl Nonlinear Dyn 2019;8(2):109–66. 17 S.R. Nekoo and A. Ollero ISA Transactions xxx (xxxx) xxx [56] Cloutier JR, D’Souza CN, Mracek CP. Nonlinear regulation and nonlinear H- infinity control via the state-dependent riccati equation technique: Part 1, theory. In: Proceedings of the International conference on nonlinear problems in aviation and aerospace. Embry Riddle University Press, Daytona Beach, Florida; 1996, p. 117–31. [57] Pang H, Wang L. Global robust optimal sliding mode control for a class of affine nonlinear systems with uncertainties based on SDRE. In: 2009 Second International workshop on computer science and engineering. 2, IEEE; 2009, p. 276–80. [58] Franzini G, Innocenti M. Nonlinear H-infinity control of relative motion in space via the state-dependent riccati equations. In: 2015 54th IEEE conference on decision and control. IEEE; 2015, p. 3409–14. [59] Feydi A, Elloumi S, Jammazi C, Benhadj Braiek N. Decentralized finite-horizon suboptimal control for nonlinear interconnected dynamic systems using SDRE approach. Trans Inst Meas Control 2019;41(11): 3264–75. [60] Nekoo SR. Output-and state-dependent riccati equation: An output feedback controller design. Aerosp Sci Technol 2022;126:107649. 18