I Introduction
The dynamic modeling of mechanical and bioinspired robotic systems is often conducted using ordinary differential equations (ODE) that describe the evolution of states such as positions and orientations over time. These models are obtained based on known physics of the mechanical systems and generally many interconnected processes inside of the described system are not known and cannot be modeled using ODE. For this reason, machine learning methods that identify the whole ODE
[ogunmolu_nonlinear_2016, lu_robust_1998] and physicsinformed machine learning methods that fill the gaps in models based on ODE framework have emerged[raissi_physicsinformed_2019, stiasny_physicsinformed_2021, qian_physics_2019]Nevertheless, these machine learning methods are limited to offline training since they are trained using backpropagation and gradient descent methods. As a result, their training methods require huge data sets to converge for a solution which limits their time performance. The time of training particularly becomes important if these models are meant to be used in conjunction with fast innerloop closedloop feedback. This particularly becomes very problematic in robotic platforms that have highly nonlinear and continuous interactions with their fluidic environment, such as aerial and aquatic locomotion. While the application of machine learning is expanding in these disciplines, yet it is too early to make conclusive decision about their suitability for realtime closedloop feedback because their accuracy and computation overhead is unknown.
In other applications related to robot locomotion, often neural networks are trained around a local attractor with no possibility of updating neural networks weights in an online fashion. To enable online training, Bayesian approaches for training neural networks have emerged [noauthor_training_nodate, wang_convergence_2011, safarinejadian_predict_2013]. However, the application of these Bayesian approaches remained limited.
In this paper, we attempt to capitalize on an ongoing research project at Northeastern University, called Aerobat [sihite_computational_2020], to explore novel methods of estimating complex robotic models using machine learning and neural net policies. Briefly speaking, Aerobat is a morphing Micro Aerial Vehicle (MAV) and its flight dynamics contains hardtomodel components involving complex fluidstructure interactions.
The objective of this paper is to boost our efforts in understanding Aerobat’s complex dynamics with a platform shown in Fig. 1. Using the platform, it is possible to define any desired trajectories for the robot manipulator and Aerobat which the former defines flight path and the latter defines robot inputs. Then, the goal is to use the robot’s Lagrangian model and measurements from onboard sensors to estimate models that realistically capture fluidstructure interactions in Aerobat.
To do this, we use Algorithmic Differentiation (AD) and, inspired by [arasaratnam_cubature_20091]
, use Kalman filters that employ cubature approximation rules to numerically compute Gaussianbased integrals. As a result, we show that we can efficiently train the physicsinformed, feedforward neural network that captures Aerobat’s dynamics.
This work is organized as follows. After briefly introducing NU’s Aerobat, we drive the equations of motion using the method of Lagrange. Next, we describe two numerical approaches based on AD and cubature rules to extract efficient numerical models suitable for modelbased nonlinear control of Aerobat. Our algorithmic computations are developed based on CasADi toolbox which is briefly described. The bulk of this works is focused on aerodynamic force estimation in our platform. We conclude this work with final remarks and simulation results.
Ii Brief Overview of Northeastern Morphing MAV
The Northeastern University’s Aerobat, which is a morphing MAV shown in Fig. 2, weighs less than 20 gr with a maximum wingspan of 32 cm and is capable of concurrently mobilizing 14 active linkage DOF in gait cycles shorter than 200 milliseconds. We have developed this platform (and the developement is still ongoing) to study control and actuation frameworks that can help us design future morphing MAVs.
Iii Morphing MAV Model
Here, we describe the equations of motion which will be used to obtain the simulation results reported in this paper. The considered flight dynamics assume interconnected inertial dynamics and aerodynamic subsystem which is expected from these systems as they do not follow insect flight dynamics. The dynamics, as will be shown below, is affineincontrol which is less common to see in flapping wing flight.
Let be the underactuated (passive) coordinates (body position and orientation parameterized with Euler angles) and
be the actuated (active) coordinate vectors (wing joint angles). These underactuated and actuated coordinates result in a system defined by the configuration variable vector
and the state vector .Next, the equations of motion are derived using the EulerLagrange formalism after obtaining the system’s Lagrangian . In obtaining the Lagrangian functional, the kinetic (translational and rotational) and potential energy led by the distributed mass from each wing segment is considered which for saving space explaining them is overlooked in this paper. Then the equations of motion are derived using the general formulation
(1) 
where is the sum of the generalized forces. This formulation can be expanded into the following form:
(2) 
where denotes the mass inertia matrix, contains the Coriolis terms, is the gravity vector, are the joint torques, maps joint torques and (similar to ) maps external forces to the generalized coordinates in the system. In Eq. 2, the aerodynamic contributions are considered as the second input to the system.
The aerodynamic forces are modeled by obtaining the resultant lift and drag forces on discrete blade elements at their quarterchord location denoted by () shown in Fig. 3. The combined aerodynamic generalized forces can be derived as follows:
(3) 
where embodies all of the quarterchord points, see Fig. 3. We model , which is a skinny vector containing all of the aerodynamic forces at each wing section, using models of form given in [izraelevitz_statespace_2017].
In this work, our objective is twopronged and include

First, we show that using algorithmic differentiation it is possible to generated models of articulated flapping flight useful for realtime modelbased flight control.

Second, we demonstrate that iterative model estimations for can be achieved in a feasible way (from a closedloopfeedbackdesign standpoint), allowing the accurate prediction of nonlinear fluidstructure interactions in realtime.
Algorithmic Differentiation of Inertial, Coriolis Effects, Gravity in Morphing Wings
In Eq. 2, the inertial matrix, , Coriolis ,, gravity vector, , are computed based on a precise modeling scheme that considers the joint angles, , their velocities, , and applies forward and velocity kinematics to obtain position and velocity terms symbolically. We obtain these terms symbolically using Matlab’s Symbolic Toolbox. Following standardized steps , and are obtained. However, these symbolic matrices are massive and their application for fast flight control and motion planning is infeasible.
The alternative to this is Algorithmic Differentiation (AD) [carius_trajectory_2018] which is as accurate as symbolically generated results but much leads to much faster computation overhead. There are two approaches to do AD: (i) operator overloading or (ii) source code transformation. We have decided to use the second approach using the toolbox called CasADi because of its great success in pathplanning and motion control in robot manipulators and floatingbase legged systems. CasADi uses a powerful parser to introduce the derivative terms.
Iv Efficient FluidStructure Interaction Model Estimation Using Cubature Rules
Here, our objective is to find estimated models for fluidstructure interactions in articulated flapping flight. The main motivation is to obtain models that have low computation overhead, can accurately predict aerodynamic contributions, and as a result can be used for realtime flight control. It is possible to show that these fluidstructure interactions, in addition to the state vector , can be parameterized in terms of some internal variables in form of indicial models given by
(4) 
where , , and are some state and timedependent matrices. The hidden variables are denoted by , and defines the holonomic constraints that describe morphing in the system (enforced by from Eq. 2). The form of Eq. 4 has motivated us to use neural network policies to find aerodynamic models described by
(5) 
where is the statedependent, internal dynamics capturing fluidstructure interactions and is a nonlinear function of the hidden variable . Note that Eq. 5 partly captures the behavior of Eq. 4 as it takes the state vector as the input and ignores ’s role in steering . So, the idea is to capture this hidden dynamics using neural net policies. In other words, throughout this section our main objective is to sample the robot’s configurations (i.e., ) at time points during flight and utilize vehicles Lagrangian dynamics, which present the predictable part of flight dynamics, to train a deep neural network such that an error norm is minimized, which will be briefly explained later.
The deep neural network is a hierarchical model of multiple layers and activation functions, where each layer applies a linear transformation followed by a nonlinear activation function
to the preceding layer. The equation of is given by(6) 
Where is the vector containing the weights of ith layer and embodies the sampled robot configurations. The activation function at ith layer is denoted by . So, the general idea is to use sampled robot configurations to construct training data and then use that to train the neural net weights such that the error norm given below is minimized
(7) 
where is the optimal neural network weight. Note that the training data (j denotes sample index) is directly obtained from the model dynamics as
(8)  
where and denote robot configurations at jth sample time. Backpropagation is usually used to find the optimal weights. Since backproppagation requires large amount of time to find , many samples (i.e., skinny matrix) have to be collected offline and the optimization process have to also be done offline.
Another problems with backprogagation involves narrow boundaries of prediciton which becomes very important in our application. Meaning to realistically predict fluidstructure interactions, the measurements must be collected at various robot configurations which is not feasible. As a result, we can end up with measurements collected and weights trained around a local attractor. There is the risk that when the robot states escape that local attractor, neural network weights cannot be updated in realtime due to the computational cost of backproppagation.
Extended Kalman Filter (EKF) can be considered to estimate [ljung_asymptotic_1979]. Since in the neural network case, EKF method depends on computing the Jacobian term the computational cost can be very high. In this paper, we use cubature numerical approximation rules to find the optimal weights . Unlike the EKF [ljung_asymptotic_1979], the cubature rules do not need to solve the jacobian term , and can offer practical numerical solutions for computing the Gaussianweighted integrals that appear in posterior terms involved in model estimation.
Iva Model Estimation Problem
In general, the KF training algorithm can be divided into two classes, parallel KF and parameterbased KF. In parallel KF, both neural network output and weights are the states to be estimated by the KF algorithm. For parameterbased KF, only weights are treated as the states to be estimated. We consider the second approach which treats the training process as a filtering process.
We will turn the problem of finding the unknown internal dynamics given by Eq. 5 into the following neural network model estimation problem. Consider the neural network weights to be found in Eq. 7. The following discritized dynamics are considered
(9) 
where the subscripts denote the kth sample time; and are the process and measurement noise, respectively. We stack all of the weight parameters and associated noise from ith layers to simplify the notation as following
(10) 
where and embody the stacked parameters from each layer. Note that the neural network weights are the state vector and the vehicle configuration is the input in Eq. 10. While, in Eq. 5, we assumed the output is a function of states and hidden variable , in Eq. 10, we assume the output function is the function of states and neural net weights () which explains how the hidden dynamics is replaced by neural net policies to be found in this section.
We also note that in Eq. 10, the main reason that the process dynamics is simply identity is because the estimation problem is solved under the assumption that the weights are optimal weights in which case their value must remain unchanged, i.e., . Ideally, must be constant. However, we add the artificial process noise to provide more flexibility in tuning the filter.
IvB Bayesian Estimation
The general form for the estimated dynamics of Eq. 10 is given by
(11) 
where is the estimated weight vector at kth time step, is the generated Kalman gain according to cubature approximation rules inspired by [arasaratnam_cubature_20091], and is the output of the neural network.
Assume we have a set of measurements , the goal of any bayesian estimator is to estimate state form measurements by finding the posterior density function
(12) 
where
denotes the probability density function and its fully determined by it is mean value and covariance matrix. The posterior density function
is fully determined by obtaining the mean and covariance matrix of , and . Since finding the mean and covariance terms for all of these terms follow a similar procedure, here, let us only show obtaining the mean value and covariance matrix of .The expected value of given the a priori information is given by the following generic equation
(13) 
where is the space on which is defined. Note that to write Eq. 13, Eq. 10 and the fact the
is zeromean random variable are used. In Eq.
13, is the Gaussian density function with the mean value given by and covariance matrix given by . The covariance matrix is given by(14)  
where is the process noise covariance matrix at kth sample time. Next, we briefly describe the approximation methods used to find Guassianweighted integrals in the above equations.
IvC Approximation of Probability Density Terms
When the predictive (process) and observation models are linear, that is the Gaussianweighted integrals to be solved possess the form
(15) 
where is a linear function of state vector and is the a Gaussian density function, then, the integral given above can be solved analytically, leading to the wellknown Kalman time and measurement update equations[calandra_bayesian_2016, rai_bayesian_2017, ho_bayesian_1964].
When the predictive and observation models are nonlinear as is the case presented by Eq. 9, the linearization of the model [ljung_asymptotic_1979] or numerical approximation solutions must be implemented. For instance, EKF linearizes the model using Taylor series expansion. The linezarization process can be computationally expensive for the application of Eq. 9 in realtime closedloop feedback flight control of our models. Mainly because, EKF heavily depends on
at every time step. To solve this issue, we consider filters that use cubature rules to compute the Guassianweighted integrals that appear in the conditional probability density function explained previously.
Consider the following nonlinear Gaussianweighted integral, which captures terms such as Eq. 14,
(16) 
where denotes any nonlinear function of . Then, the integral given by Eq. 16, can be approximated by
(17) 
where , and . Let (where and are …) and , then it is possible to show that
(18) 
(19) 
which is directly used to write the algorithm given below to compute the approximation of Gaussianweighted integrals in the Kalman Filter.
V Simulation Results and Discussion
Here, we briefly report the accuracy and computation overhead of the obtained model for our morphing wing flight based on AD and cubature rules. Our results suggest the feasibility of using these models for locomotion control of robots whose dynamics involves sophisticated robotenvironment interactions.
Figures 4 and 5 show the computation overhead for various terms in the inertial dynamics. As expected AD yields significantly smaller computation overhead. To check the accuracy of AD models versue symbolically generated terms, we use ADbased inertia, Coriollis and gravity terms to obtain training data points for our cubature filter.
Note that we use a mathematical model to extract our training data . However, the proposed framework is developed with the aim of implementing it on to the actual robot. So, without loss of generality the same training data can be obtained by capturing the vehicle’s body position, orientation, wing angles and their times derivatives and using the vehicles Lagrangian dynamics (excluding aerodynamics interactions).
Figures 6 and 7 show a comparison between the estimated and actual models. In these simulation results, our neural network consists of three layers with two softplus activation functions for form . At every 100 samples, one training data is utilized towards obtaining optimal neural network weights using the cubature Kalman Filter. And then the learned models is used to estimated fluid structure interactions in the robot.
The robots’ body joint angles are commanded with predefined values. Damping coefficients are considered to passively stabilize roll, pitch, and yaw dynamics as the vehicle is openloop unstable and no controller is involved.
Based on our simulation results, and as it can be seen in Figs. 6 and 7 the estimated aerodynamics model can do a decent job predicting actual generalized forces in a 200milliseconds time envelope (sample time is sec). The process and measurement noise values ( and ) are assumed to be zero. The mean square error is for all generalized forces. The diagonal terms in the covariance matrix are less than which reflects the confidence of the filter about the estimated weights.
Vi Concluding Remarks
In this work, algorithmic differentiation and neural networks equipped with cubature rules were proposed to identify inertial and fluidic forces acting on a morphing MAV called Aerobat. These forces can be very complex to the extent that averaged models fail to predict them. The neural network design extends to prior works on Bayesian filters. Using cubature rules to compute Gaussianweighted integrals efficiently, we showed that the complex multidegreesoffreedom dynamics of morphing MAVs can be computed very efficiently and accurately. In the future, we aim to integrate online fluidic structures identification of these morphing MAV with realtime closedloop feedback control and test them experimentally on our Aerobat platform.
Comments
There are no comments yet.