The dynamic modeling of mechanical and bio-inspired 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 physics-informed machine learning methods that fill the gaps in models based on ODE framework have emerged[raissi_physics-informed_2019, stiasny_physics-informed_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 inner-loop closed-loop 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 real-time closed-loop 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 hard-to-model components involving complex fluid-structure 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 fluid-structure interactions in Aerobat.
To do this, we use Algorithmic Differentiation (AD) and, inspired by [arasaratnam_cubature_2009-1]
, use Kalman filters that employ cubature approximation rules to numerically compute Gaussian-based integrals. As a result, we show that we can efficiently train the physics-informed, feed-forward 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 model-based 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 affine-in-control 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 vectorand the state vector .
Next, the equations of motion are derived using the Euler-Lagrange 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
where is the sum of the generalized forces. This formulation can be expanded into the following form:
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 quarter-chord location denoted by () shown in Fig. 3. The combined aerodynamic generalized forces can be derived as follows:
where embodies all of the quarter-chord 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_state-space_2017].
In this work, our objective is two-pronged and include
First, we show that using algorithmic differentiation it is possible to generated models of articulated flapping flight useful for real-time model-based flight control.
Second, we demonstrate that iterative model estimations for can be achieved in a feasible way (from a closed-loop-feedback-design standpoint), allowing the accurate prediction of nonlinear fluid-structure interactions in real-time.
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 path-planning and motion control in robot manipulators and floating-base legged systems. CasADi uses a powerful parser to introduce the derivative terms.
Iv Efficient Fluid-Structure Interaction Model Estimation Using Cubature Rules
Here, our objective is to find estimated models for fluid-structure 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 real-time flight control. It is possible to show that these fluid-structure interactions, in addition to the state vector , can be parameterized in terms of some internal variables in form of indicial models given by
where , , and are some state- and time-dependent 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
where is the state-dependent, internal dynamics capturing fluid-structure 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.
Where is the vector containing the weights of i-th layer and embodies the sampled robot configurations. The activation function at i-th 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
where is the optimal neural network weight. Note that the training data (j denotes sample index) is directly obtained from the model dynamics as
where and denote robot configurations at j-th 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 fluid-structure 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 real-time 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 Gaussian-weighted integrals that appear in posterior terms involved in model estimation.
Iv-a Model Estimation Problem
In general, the KF training algorithm can be divided into two classes, parallel KF and parameter-based KF. In parallel KF, both neural network output and weights are the states to be estimated by the KF algorithm. For parameter-based 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
where the subscripts denote the k-th sample time; and are the process and measurement noise, respectively. We stack all of the weight parameters and associated noise from i-th layers to simplify the notation as following
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.
Iv-B Bayesian Estimation
The general form for the estimated dynamics of Eq. 10 is given by
where is the estimated weight vector at k-th time step, is the generated Kalman gain according to cubature approximation rules inspired by [arasaratnam_cubature_2009-1], 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
denotes the probability density function and its fully determined by it is mean value and covariance matrix. The posterior density functionis 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
is zero-mean 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
where is the process noise covariance matrix at k-th sample time. Next, we briefly describe the approximation methods used to find Guassian-weighted integrals in the above equations.
Iv-C Approximation of Probability Density Terms
When the predictive (process) and observation models are linear, that is the Gaussian-weighted integrals to be solved possess the form
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 well-known 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 real-time closed-loop 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 Guassian-weighted integrals that appear in the conditional probability density function explained previously.
Consider the following nonlinear Gaussian-weighted integral, which captures terms such as Eq. 14,
where denotes any nonlinear function of . Then, the integral given by Eq. 16, can be approximated by
where , and . Let (where and are …) and , then it is possible to show that
which is directly used to write the algorithm given below to compute the approximation of Gaussian-weighted 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 robot-environment 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 AD-based 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 open-loop 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 200-milliseconds 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 Gaussian-weighted integrals efficiently, we showed that the complex multi-degrees-of-freedom 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 real-time closed-loop feedback control and test them experimentally on our Aerobat platform.