Annual Reviews in Control Decoding and realising flapping flight with port-Hamiltonian system theory

In this paper we envision how to tackle a particularly challenging problem which presents highly interdisci- plinary features, ranging from biology to engineering: the dynamic description and technological realisation of flapping flight. This document explains why, in order to gain new insights into this topic, we chose to employ port-Hamiltonian theory. We discuss how the physically unifying character of the framework is able to describe flapping dynamics in all its important aspects. The technological and theoretical challenges of flapping flight are discussed by considering the interplay between different topics. First of all, the formal conceptualisation of the problem is analysed. Second, the features and capabilities of port-Hamiltonian framework as the underneath mathematical language are presented. Subsequently, the discretisation of the resulting model by means of structure-preserving strategies is addressed. Once a reliable numerical model is available, we discuss how control actions can be computed based on high-level specifications aiming at increasing the flight performances. In the last part, the technological tools needed to validate experimentally the models and to equip a robotic bird prototype with the necessary sensing and actuation devices are discussed.


Introduction
Flapping flight is one of the wonders of nature and has been studied vastly by biologists and fluid dynamicists. The level of dexterity, efficiency and grace that birds exhibit along flight manoeuvres is far beyond what state-of-the-art technology can achieve. Furthermore, a complete description able to disentangle the complexity of the emerging fluid dynamic patterns and to yield an interpretation of all flight modes, is not covered by current theoretical understanding. Even if this lack of knowledge is recognised worldwide, it remains astonishing, especially if one thinks that the basic principles which play a role in flight are rooted in classical physics of continua. The reason is that classical physics can be very mysterious when it comes to fluid dynamics, especially in those regimes where viscosity, turbulence and vortex generation are important. These are exactly the aerodynamic phenomena that birds use in flapping flight, generating lift and thrust by means of sophisticated generation of unsteady aerodynamic patterns interacting on their wings (Chin & Lentink, 2016).
In physical terms, and supposing that no principle is missing in our description of continua, understanding bird flight would follow from understanding the time evolution of an extremely complicated fluid-solid interaction system: on one side the fluid (whose evolution is described by 3-D Navier-Stokes equations), and on the other the wing (described by 3-D elasticity equations) play their role in an entangled the port-Hamiltonian methodology forward, it is possible to handle dynamic boundary conditions as fundamental part of the system, while preserving the true power flow between the two continua in the game, without simplifying the underlying physical models. In this way we expect to gain new insight in fluid patterns and their evolution along unsteady solutions.
The port-Hamiltonian structure of the model can be used to build numerical integration schemes able to preserve the important physical properties of the underlying continuous system, with the consequent possibility of generating physical reliable data in dynamic conditions relevant to flapping flight. We show how the model can then be validated by wind tunnel tests with flow visualisation, and how numerical optimisation will play a role in tuning control inputs which will realise different flight modes. Finally, after briefly reviewing the current approaches used to design robotic artefacts, we envision how a new robotic bird can be realised with unprecedented flight dexterity, able to flap asymmetrically, adapt to the flow and take off and land as birds do, in order to validate the scientific understandings.
This paper is organised as follows. In Section 2 an overview of the important mechanisms involving flapping flight regimes is given, which will serve as motivation to claim the necessity of a theoretical paradigm shift with respect to standard approaches. In Section 3 the core methodology -port-Hamiltonian theory -that will be used to model the system, is presented. Section 4 contains a description on how novel numeric methods can be developed by taking advantage of the developed model. Section 5 discusses how the control problem will be tackled and Section 6 discusses possible technological solutions and mechanical design principles that will be implemented to finally realise a prototype robotic bird embodying the discussed principles.

Modelling flapping flight
The understanding of the aerodynamic, structural, and behavioural patterns of birds is essential for the development of robust and high performance flapping-wing aerial robots. However, this is a tremendously difficult task.
In order to grasp the complexity that birds deal with, we give some examples. They (i) exploit multiple aerodynamic mechanisms for the generation and enhancement of lift and thrust; (ii) accommodate wind gusts and avoid obstacles via variable kinematic and dynamic patterns using their wings, tail, and body; (iii) actively change the structural compliance of their wings to achieve wing deformations that improve aerodynamic performance; and (iv) have distributed actuation and sensing biological units all over their surface that are used to ensure flight stability in extreme environmental regimes (Shyy, Aono, Kang, & Liu, 2013).
Despite the many research efforts throughout the past decades, flapping flight dynamics is still not fully understood (Floreano, Zufferey, Srinivasan, & Ellington, 2009;Shyy et al., 2013), especially that of bats and birds (Chin & Lentink, 2016). The latter reference, in particular, presents an excellent state-of-the art survey about flight mechanisms ranging from insects to vertebrates. It is evident that even if a number of known phenomena can be disentangled to understand aerial locomotion, e.g. leading edge vortices (LEV) to prevent stall, we are far from having a physical description of the mechanism underlying the important phenomena. For example, there is not a general 3D modelling framework which is able to reproduce the generation of the LEVs in all situations of both sustained motion, take-off and hovering and can give conclusive results and explanations of what is observed experimentally.
What is evident is that to properly describe in a model the flappingflight of birds, essential ingredients that need to be taken into account without simplifying assumptions are viscosity of the airflow and the generation of vortices. These features are properly described in an unsteady 3D fluid-solid interaction system in which bilateral wake-structure interactions, adaptive wing morphologies, and highly nonlinear fluid dynamic effects play a significant role. From a dynamic modelling perspective, unless considering degenerate cases, a fluid-solid interaction system is not fully determined by the dynamic evolution of the two parts separately: The so called interface conditions transform the interconnected parts into a single, completely new dynamic system which can exhibit rich and chaotic behaviour at the interface. As consequence, escaping from this complexity using patches of quasi-steady aeroelastic models means surrendering in trying to decode flapping flight principles. In fact we know that flight exploits wake patterns that cause dynamic interaction back to the wing, a phenomena not captured in quasi steady models, but in the full 3D Navier-Stokes equations. Unsteady aerodynamic models able to predict nonlinear effects (e.g., flutter, edge suction, etc.) exist and have been demonstrated to be very useful for analysis and validation against numerical methods (see e.g. Boutet and Dimitriadis (2018) for fixed wing and DeLaurier (1993) for flapping wing). Nevertheless their validity is limited to specific ranges of kinematic and dynamic parameters of the system, and it is extremely difficult to make use of these models to describe complex motions with large kinematic variations (e.g. high angles of attack) and to design model-based control laws, since the controller could structurally modify the validity range of the models themselves.
In any case, the fact that birds fly indicates that, even if the complexity of this system is very high and the fluid exhibits turbulence (and thus the whole system possesses a chaotic evolution) it is possible to learn dynamic patterns and exploit them.
We accept this challenge using a new complementary framework which could give new insight in fluid patterns and their evolution. The claim is that a physically reliable model of the fluid-solid system can reveal truth about flapping flight mechanism, overcoming the difficulty in dealing with the chaotic nature governing the dynamics. By physical reliability we mean that the model must be able to keep track, in a coordinate-free way, of the energy exchanges within the system while respecting the real physical conservation laws characterising the modelled processes. As an example, this change of perspective was successfully applied in rigid-body dynamics, giving rise to the recent energy-aware robotics framework (Stramigioli, 2015), which led to novel control strategies, outlined in Section 5.
In this context the methodology needs to be pushed forward in order to deal with coupled, spatially distributed systems. To summarise, we search for a modelling paradigm which is able to keep track of the energy exchanges within the system parts preserving their nonlinear nature, and that systematically takes into account dynamic boundary conditions in the fluid-solid interaction system. In the next section we discuss why the port-Hamiltonian (pH) paradigm is perfectly designed to achieve this goal.

The port-Hamiltonian framework
In this section we start by briefly giving a ''bird's eye'' view on the main features characterising the paradigm of port-Hamiltonian theory underlying physical systems, and then apply them to the case of modelling a Bio-inspired robotic bird.

The Geometric port-Hamiltonian approach
The port-Hamiltonian (pH) framework (Maschke & van der Schaft, 1992) aims at modelling and controlling physical dynamical systems in a substantially different way than the majority of the control system community does. It uses concepts like power flows and energy to constrain the model to follow physical laws, rather than of signals, which can be seen as the quantities ''flowing along the arrows'' of conventional block diagram representation of I/O dynamical systems. It is rooted in two classical frameworks: 1. The port-based modelling paradigm, introduced by H. Paynter in the Bond-Graph formalism (Filippo, Delgado, Brie, & Paynter, 1991), which describes a complex system as the interconnection of simpler atomic units exchanging energy.
F. Califano et al. Fig. 1. A port-Hamiltonian system consisting of an energy storage subsystem, an energy dissipation subsystem, an external interaction port, in addition to a Dirac structure that acts as an energy routing element.
2. The Hamiltonian mathematical formalism (Marsden & Ratiu, 2013), which is one of the two main methodologies (the other is the Lagrangian approach) to model physical systems, and it emphasises the geometric role of the phase space and of the Hamiltonian function (representing the total stored energy) as building blocks for describing mechanical systems. In practice, using this formulation only closed, or isolated systems are modelled, i.e. systems for which the Hamiltonian is a conserved quantity along the motion.
The pH framework elegantly synthesises the two ideas: it extends the geometric Hamiltonian formulation to open systems of any physical domain, i.e. to systems that follow classical Hamiltonian dynamics if isolated (if the input is set to zero), but are equipped with collocated input-output pairs. These allow the interconnection with other systems and the possibility to incorporate energy dissipation in the framework. It becomes thus possible to model increasingly large networks in which energy transfers between subsystem are systematically monitored and possibly used for control purposes. This energetic modular property is not shared with classical I/O block diagram modelling approaches which, contrarily to pH systems, (i) do not reflect the physical nature of conservation laws (do not present an autonomous Hamiltonian structure), (ii) present a fixed and arbitrarily causality choice from the beginning (there is a fixed and non physical choice of input and output). As comparison pH models are acausal, in the sense that each subsystem contains only the constitutive relations describing how its variables evolve without imposing which is an input and which is an output. The single constitutive relations of the system are classified on the basis of their energetic behaviour into energy-storing, energy-dissipating, energysupplier, or energy-routing systems. The latter type encodes the power continuous network structure of the system, describing how the power flows are constrained, and is mathematically modelled by the so called Dirac structure, generalising the Poisson structure, peculiar of classical Hamiltonian treatment.
The presence of nonlinearities in the system is encoded either in the constitutive relations or in state-modulated Dirac Structures. This is one of the most powerful properties of pH systems: there is no conceptual increase of difficulty when nonlinear systems need to be modelled. This constitutes a major difference with respect to classical control theoretic approaches, which tend to classify the nature of the systems on the basis of the form of the equations, rather than on the underlying geometric properties.
An interaction between two systems is characterised by the effect of the energy exchanged between them. This interaction takes place through the aforementioned collocated I/O pairs, which constitute a power port. Each power port consists of two dual variables, called an effort and flow, whose duality pairing represents the power flowing between the two interacting systems. We refer to Fig. 1 to aid the reader with a graphical representation, using the Bond-Graph language, of the discussed concepts.
By defining these ports as natural, geometrical, duality structures which are present in all domains in physics and using them as the basis for interconnection, it possible to naturally describe systems of any kind and dimension as the natural interconnection of parts or phenomena.
The mathematics in which the framework is developed uses the coordinate-free language of differential geometry rather than vector calculus, which depends on a particular coordinate system. This provides an elegant generalisation of pH theory to smooth manifolds and to infinite-dimensional systems (Van Der Schaft & Maschke, 2002) (see Rashad, Califano, van der Schaft, and Stramigioli (2020) for a recent survey). In fact pH systems are structurally suited to handle distributed systems with boundaries through which non-zero energy is exchanged. This is a fundamental feature in the modelling of flapping, due to the interaction between the flexible mechanical domain and the aerodynamic domain. Furthermore, in contrast to classical approaches, the framework can be used to represent interconnections of finite and infinite-dimensional systems. Such capability has been demonstrated e.g. in modelling robotic manipulators with mixed rigid and flexible links (Macchelli, Melchiorri, & Stramigioli, 2009) and interconnected distributed and lumped systems (e.g. a flexible plate welded to a rigid bar (Brugnoli, Alazard, Pommier-Budinger, & Matignon, 2019)).

Port-Hamiltonian view of a Bio-inspired Robotic Bird
In the port-Hamiltonian framework, a flapping bird, viewed as a dynamical system, can be abstractly modelled as a number of interconnected subsystems, as shown in Fig. 2.
First, the main hull of the bird can in first approximation be modelled as a single rigid body interconnected to two wings via joints and power ports. Second, the two wings are viewed as two flexible bodies with distributed stiffness and mass corresponding to the different components of the biological wing such as bones, muscles, and feathers. Along the wings, distributed actuators model the mechanisms that birds utilise to actively deform their wings, and distributed sensors model the mechanisms that birds use to, for example, sense air flow. Third, both the bird's rigid body and flexible wing models are interconnected to another subsystem that represents the dynamics of the air flow. This latter interconnection, implemented by means of power ports, constitutes a fundamentally new way to construct the model of a fluidsolid interaction system. In the case of a robotic bird, the overall port-Hamiltonian model would include additional subsystems such as the control system interacting, via an electrical port, with the robot's flapping mechanism and distributed actuators engineered to deform the wings or change its stiffness. As an example, the left part of Fig. 3 gives the port-Hamiltonian view of the joint connecting the main hull and the wing of the bird.
As outlined at the beginning of this section, the port-Hamiltonian framework has a number of significant features that make it possible to model and decompose bird's flapping flight into the set of smaller subsystems described above. In particular the subsystems shown in Fig. 2 are a combination of finite and infinite-dimensional systems: the rigid body model is a finite-dimensional system represented by ordinary differential equations, whereas the airflow and flexible body models are infinite-dimensional systems represented by partial differential equations.
Furthermore, as described in Section 2, to properly describe the mechanism of flapping-flight, essential aerodynamic phenomena that need to be considered are viscosity of the air and the generation of vortices. Consequently, linear or ideal assumptions cannot be made in the mathematical modelling of the air and the full Navier-Stokes equations have to be used, which makes the airflow model a nonlinear and dissipative dynamical system. Although the traditional Hamiltonian framework, based on Poisson structures, can handle nonlinear infinite-dimensional systems, it is limited to conservative systems. The pH framework overcomes this limitation by exploiting dissipative ports that allows handling non conservative phenomena, i.e. systems that dissipate free energy. Therefore, the Navier-Stokes equations can be described fully in the port-Hamiltonian framework by interconnecting a   power preserving fluid dynamical system (referred to as Euler equations in the PDE context) to a dissipative system encoding the constitutive relations of stress in Newtonian fluids. In this spirit, lets ''zoom in'' the fluid dynamic model to conceptually show how this modular interconnection gives rise to a new way to represent Navier-Stokes equations. The right part of Fig. 3 gives an intuition about the description of the pH model for viscous fluid systems as an interconnection of a conservative part (see Stramigioli (2021a, 2021b) for early results) and a dissipative unit encoding constitutive relations for Newtonian stress (Califano, Rashad, Schuller, & Stramigioli, 2021). Even if the mathematical details to make the interconnection precise are not trivial, conceptually and from a bondgraph perspective, the procedure is straightforward. To give a grasp, it works in the same way as the interconnection of an electrical resistance within a lumped LC circuit, i.e. adding an energy dissipating element to a conservative system. Contrarily to what happens if one works with equations only, this methodology is convenient if different/additional constitutive relations need to be used, since the network structure remains unchanged and one can modify the subsystem of interest only, providing a truly open approach for network modelling. Furthermore, using the language of differential geometry, Navier-Stokes equations can be elegantly generalised on a (possibly curved) manifold, together with their energetic properties.
In order to describe a fluid dynamic system or a flexible structural system as open systems, which are able to be interconnected to other systems, the variable (spatial) boundary conditions need to be incorporated. This is the case for a fluid-solid interconnection problem, or a multi-fluid system, in which the dynamic interaction between the systems in the game induces a change in time of their respective spatial domains. This aspect is often referred to as moving interface, and a proper geometric treatment of the problem in the pH framework is under current study.
A major limitation of the traditional Hamiltonian treatment of such infinite-dimensional systems is that it handles only boundary conditions that cause zero-power exchange through the boundary. Thus, it can only represent isolated systems which are not part of a bigger dynamical process. In contrast, as highlighted previously, the pH paradigm can properly handle dynamically changing boundary conditions by defining a power preserving interconnection between the two spatially adjacent systems. This constitutes a major generalisation with respect to PDEbased approaches, in which boundary conditions need to be defined separately and are very difficult to handle when a dynamic relation constraints these conditions, e.g. in a fluid-solid interaction problem.
To conclude, the port-Hamiltonian framework is capable of incorporating multi-domain physical systems since it is based on energyconservation principles. As a result, the same mathematical representations and concepts can be used for the different physical domains available in the model of flapping-flight i.e. structural mechanics, and fluid mechanics. In addition, such a property is also useful for modelling of distributed sensors and actuators that will lie in the electromagnetic domain. This has been demonstrated previously in modelling a complex multi-physical system such as thermo-magnetohydrodynamics (Vu, Lefevre, & Maschke, 2016) with coupled fluid mechanical, electromagnetic, and thermodynamic components.
We believe that the pH framework provides a complementary way with respect to classical methods for studying the system at hand. This modelling methodology deeply relates an engineering approach based on interconnecting simple systems, to the physical truth of prediction that the model would provide, due to its geometric and energy-based nature. In this sense, and with the help of discretisation methods described in the next section and experimental model validation, it is reasonable to believe in gaining new insights on aerial locomotion.

Structure-preserving discretisation
The considered problem is a very challenging one under a computational viewpoint. Flapping flight is a completely unsteady, multiphysics phenomenon, where the fluid-structure interaction takes place on a moving boundary. The discretisation method aims to be structurepreserving (Christiansen, Munthe-Kaas, & Owren, 2011), meaning that the main features of the distributed system (e.g. passivity or energyconservation, symmetries, reversibility, collocated input-output behaviour) are retained at the discrete level. The employment of a structure-preserving discretisation strategy is motivated by the fact that symplectic (i.e. structure preserving) algorithms do not only preserve geometrical structures of the flow of the differential equation, but they also have favourable properties concerning their global error when the integration is performed over a very long time (Calvo & Hairer, 1995). In more quantitative terms, the global error grows linearly for symplectic numerical schemes based on the Hamiltonian model, whereas it grows quadratically for one-step methods (e.g. Runge-Kutta methods) (Hairer, 1997). Hamiltonian integrators are the method of choice for orbital mechanics, molecular dynamics and plasma physics. Our expectation is that the Hamiltonian description of flapping flight will lead to high-fidelity numerical schemes even in this complex multiphysics application. Additionally, the discretisation may take advantage from the modularity of pHs, so that each physical subsystem can be discretised separately from the others, avoiding the need for a monolithic resolution of the problem. The power preserving interconnection implemented at a continuous level in the pH model can thus lead to numeric schemes able to structurally encode the physical power flow between subsystems at a discrete level, which is a property that can hardly be achieved using standard interpolation methods to couple FEM methods for structural dynamics and CFD methods for fluid dynamics. The discretisation has to account for the semi-discretisation in space and for time integration of the resulting ODE. These two topics are discussed separately in the next sections.

Semi-discretisation in space
Structure-preserving semi-discretisation of distributed pHs is an active research topic since almost 20 years. The first study dates back to Golo, Talasila, Van Der Schaft, and Maschke (2004), where the authors made use of a mixed finite element spatial discretisation for 1D and 2D hyperbolic system of conservation laws. A pseudo-spectral method relying on higher-order global polynomial approximations was studied in Moulla, Lefevre, and Maschke (2012). This method was used and extended in Cardoso-Ribeiro, Matignon, and Pommier-Budinger (2017) to take into account unbounded control operators. A simplicial discretisation based on discrete exterior calculus was proposed in Seslija, van der Schaft, and Scherpen (2012). A 2D finite difference method with staggered grids was used in Trenchant, Ramírez, Le Gorrec, and Kotyczka (2018). Weak formulations based on finite elements began to be explored in the last years. In Kotyczka, Maschke, and Lefèvre (2018) the prototypical example of hyperbolic systems of two conservation laws was discretised by a weak formulation. Therein the boundary is split according to the causality of boundary ports, so that mixed boundary conditions are easily handled.
Recently, it has become evident that there is a strict link between discretisation of port-Hamiltonian (pH) systems and mixed finite elements (Cardoso-Ribeiro, Matignon, & Lefèvre, 2020). Velocity-stress formulation for the wave dynamics and elastodynamics problems are indeed Hamiltonian and their mixed discretisation preserves such a structure (Kirby & Kieu, 2015). Mixed finite elements allows achieving a structure-preserving discretisation even in non-linear examples, as demonstrated in Bauer and Cotter (2018) for the non-linear rotational shallow water equations. However, some form of upwinding is required to guarantee the stability of the scheme in this case.
A mixed finite element strategy represents a promising solution for the problem under consideration. There exists a strong interplay between (mixed) finite elements, differential complexes and differential geometry (Arnold, Falk, & Winther, 2006). Indeed, ''the stability and convergence of numerical methods rely on the preservation of the underlying structures of the differential complexes'' (Arnold & Hu, 2020).

Time integration of the resulting ODE
As previously stated, geometric integration of ODEs is particularly suited to perform long-time accurate simulations. A fundamental reference book on this topic is the one by Lubich and Hairer (Hairer, Lubich, & Wanner, 2006). Therein, only closed Hamiltonian systems that do not interact with an external environment are accurately studied. Given the multiphysics nature of flapping flight, it is crucial that the energy exchanges between subsystems are correctly handled by the discretisation. By using collocation methods and their associated quadrature formulas it is possible to consistently approximate both the solution and the supplied energy in the case of open Hamiltonian systems (Kotyczka & Lefèvre, 2019). For this reason, collocation methods represent a promising strategy that may be capable of handling this complex multiphysics problem in a modular manner.

Control in the pH framework
In addition to modelling, the port-Hamiltonian framework comes together with peculiar techniques for analysing and controlling complex physical systems. The interconnection structure and the total energy (i.e. Hamiltonian) of the pH system, are directly used for the design of controllers (Duindam, Macchelli, Stramigioli, & Bruyninckx, 2009). The control paradigm in the pH framework converges naturally to passivity based control (PBC) methods, which flourished under different forms, based on energy shaping and damping injection. The main idea behind these methods is to use insights on the information contained in the pH model of the system (e.g., explicit energy balance, interconnection structure) in order to (i) design control actions producing a desired closed-loop system with some different energetic properties than the uncontrolled one (energy shaping), and (ii) do it in a stable way (damping injection) (Ortega, Van Der Schaft, Mareels, & Maschke, 2001). These ideas allow in principle the analysis and control of both linear and nonlinear, lumped and distributed systems using the same fundamental concepts, and generalising the methodology to bigger class of systems, e.g. non collocated ones.
Another distinguishing feature of the port-Hamiltonian control approach is that the control system can be designed as another virtual physical system interconnected to the actual system via power ports (Stramigioli, 2001), in a methodology called control by interconnection (CbI). The advantages of such a physical interpretation of a controller are numerous. In fact the interconnection of port-Hamiltonian systems is again a port-Hamiltonian system (Cervera, van der Schaft, & Baños, 2007). Such a powerful result allows assessing the properties, like e.g. passivity of the closed loop system, by inheriting properties of its subsystems. Second, such interpretation may suggest simple and robust solutions for controlling energy flows between the system and its surrounding environment, which led to novel concepts such as a port-based interpretation of impedance control (Hogan, 1985), energy routing and energy tanks (Duindam & Stramigioli, 2004). Since the controller in the CbI paradigm resembles a virtual but physical system (e.g. a virtual spring), the design provides a physical interpretation of the control action, which could be realised by modifying the physical properties of the actual system. An example is given by the energy aware design of variable stiffness mechanisms, modulated in such a way to achieve optimal energy transfers able to steer the system on the desired behaviour. It is reasonable to believe that such a design methodology will offer new perspectives on realising optimal periodic control actions for a bird performing flapping flight, also confirmed by the biological evidence that birds do modify their stiffness parameters along a flapping period in order to maximise efficiency. Such physical implementation of the designed control action, whose feasibility makes of course contact with technological challenges addressed in Section 6, could reduce the number of sensors and actuators and lead to more energy-efficient systems, which is a fundamental requirement in the field of aerial robotics.

Learning methods from physically reliable data
Besides the previous considerations, it is a fact that the control framework developed in the context of infinite-dimensional pH systems is still too young to claim the possibility of achieving a satisfactory control design for a complex system like a flying bird using only existing model-based techniques. This is partially due to the fact that any actuation strategy would, from a control theoretic perspective, resemble an extremely under-actuated controlled system, for which most of the developed techniques are unsuitable. Furthermore on the fluid dynamic side, the vortex generation process, which is relevant to flight, is not connected to a clear sensing strategy. Even worse, also on the purely theoretical side, it is not simple to attack the control problem because the physical mechanisms determining flight modes are not precisely known and as such there is no clear link between the high level control specifications and the underlying dynamics. This is the main reason why trial and error approaches have been preferred in the past as the driving method for developing mechanisms exploiting unsteady aerodynamic effects in flapping devices. In other words, even in the presence of a physically reliable model of the fluid-solid interaction system emerging in bird flight, to take advantage of it for control purposes is not a trivial task because of both ignorance of the biological counterpart and technological limitations. In the following we explain why we believe in the heavy theoretical investment to build a physical pH model of the system, besides of the pure satisfaction of achieving a model of fluid-solid interaction with unprecedented physical reliability.
Given a fluid-solid interaction model, the control problems of interest can be cast into hard optimisation problems. For example it is straightforward to write down a problem in which we aim at maximising some high-level performance index along a flapping period (e.g. lift and thrust generation over a flapping period) while bounding the power injection from the actuators to the system (energy efficiency). These problems are nevertheless very complex to be solved because of the high dimensionality of the fluid, nonlinearity of the NS equations, and consequent non convexity of the optimisation. We stress that this difficulty in solving in a model-based fashion optimisation problems involving flapping flight in viscous fluids is ultimately the main reason why engineering solutions involving fluid dynamics historically prefer to adopt a trial and error methodology with the aid of simplified aerodynamic models. The exciting evidence, coming mainly from results developed in the last five years, is that machine learning methods are becoming exponentially better in solving hard optimisations problems using data. Even better, optimisation problems involving fluid dynamics, are embodying the perfect physical example of this synergy between hard optimisation problems and increasing availability of data coming from experiments and CFD simulations. This match, even if very recent, is so promising that a recent survey (Brunton, Noack, & Koumoutsakos, 2020) is already becoming a cornerstone in the framework.
As explained in Section 4, the pH model in its geometric formulation can be used to construct novel numeric integration schemes which, once cast in a proper simulation environment, can generate potentially unbounded amounts of physically reliable data coming from the fluidstructure interaction problem. As a consequence, it is reasonable to believe that learning methods not treated in a black box way, but using a grey box approach, i.e. the use of neural networks which contain prior knowledge of the fluid flow symmetries (e.g. conservation of energy, which is directly at the core of the pH model), can be used to gain insights on model reduction strategies, at the service of dominant pattern extraction. This idea has been successfully applied to learn a turbulence model while respecting physical properties of the fluid flow in e.g. Ling, Kurzawski, and Templeton (2016), a very recent but already important result in the field. In this respect we also refer to Duraisamy, Iaccarino, and Xiao (2019), a survey about methods to tackle closure modelling problems for turbulence models using data.
The same idea is at the core of the very recent and self explanatory Physically informed neural networks (PINN) designs on PDEs (Raissi, Perdikaris, & Karniadakis, 2019), where the neural networks are used as functional approximators of some user defined degree of freedom in the system, without destroying the physical nature of the system along the learning phase. Although admittedly challenging, there is no theoretical limitation in pushing forward these ideas in order to learn optimal actions solving the control problem at hand. In fact, a proper parametric inclusion in the model of ''physically informed'' control inputs (e.g. actuated control surfaces or variable stiffness mechanisms) would finally lead to a rational and physical-based design choice of the actuation and design principles discussed in the next section.

Creating Robots to prove Theory
A lot of theories and hypotheses can be generated, but the only real proof of having understood a certain theory and phenomenon in biomimetic locomotion is being able to synthesise it and reproduce it by creating robots and physical models which display the same dynamics. In this respect robots and physical models can even be used to set new hypotheses for biological research, especially for locomotion systems where interaction with the environment is important, like flapping flight. This approach is called robotics-inspired-biology (Gravish & Lauder, 2018). The leading edge vortex in flapping flight is an example that was first studied and understood through physical models (Dickinson, Lehmann, & Sane, 1999;Ellington, Van Berg, Willmott, & Thomas, 1996), before it was identified as a key feature for biological flight in many flying organisms (Gravish & Lauder, 2018).
There are multiple other examples of successful studies involving biological flight that have been demonstrated thanks to implementations in robotics, as shown in Fig. 4. Among the notable results we cite the following: untethered flight while matching the thrust efficiency of similarly sized insects with the Robobee (Jafferis et al., 2019); in-flight rapid and robust, underactuated wing morphing with real feathers to replicate bird wings in the Pidgeonbot ; in-phase asymmetric flapping by the Robo Raven (Gerdes et al., 2014); rapid banked turns similar to maneuvers of flies in the Delfly Nimble (Karásek et al., 2018), mimicking bat flight with articulated, flexible membrane wings in the Bat Bot (Ramezani et al., 2017); mimicking avian use of wing and tail morphing for achieving both aggressive flight and fast cruising in the LisHawk drone (Ajanic et al., 2020); and implementing active wing torsion by Festo's SmartBird (Send et al., 2012). The Robird in Fig. 4 is another example, matching its flapping dynamics and aerodynamics to that of the peregrine falcon in order to achieve similar performance (Folkertsma et al., 2017). By means of two spars and proper wing flexibility it achieves the required flapping motion. It is the first ornithopter with matching flapping dynamics and aerodynamics in the same weight and speed range of actual predator birds.
As a matter of fact, these robots produced in the last years, are moving the biomimetic robotic research field from a small set of niche engineering applications to a well established scientific branch in robotics and mechatronics (de Croon, 2020;Han, Hui, Tian, & Chen, 2020;Lau, 2020).
Despite the advancements of the state-of-the-art, every example entails only one or few solutions (e.g. wing and tail morphing, asymmetric flapping, active wing torsion, outdoor flight capabilities). The absence of a complete theoretical framework built from first physical principles advantaged a trial and error methodology to design brilliant proof-ofconcept devices but negatively affects the development of a ''general purpose'' flapping robot. In this sense, booth scientific challenges on modelling the nonlinear fluid-solid interaction system as well as engineering challenges on manufacturing proper hardware still have to be tackled in order to mature flapping wing drone technology (de Croon, 2020). We discussed how port-Hamiltonian theory would represent a theoretical framework providing a general description of bird flight in multiple regimes. In Section 5 we discussed how the methodology is  (Ramezani, Chung, & Hutchinson, 2017), e. the LisHawk (Ajanic, Feroskhan, Mintchev, Noca, & Floreano, 2020) Credit: 2020 EPFL/Alain Herzog, f. the Pidegeonbot (Chang, Matloff, Stowers, & Lentink, 2020) Credit: Lentink Lab/Stanford University, g. the Robird (Folkertsma, Straatman, Nijenhuis, Venner, & Stramigioli, 2017) and h. the Robo Raven (Gerdes et al., 2014).

Fig. 5.
In order to influence the wing-air interaction, stiffness modulation and flow sensing will be embedded. consistent with a rational and systematic choice of the sensing and actuation strategies. In the following we describe more in detail the most interesting instances of these technological solutions and then we give an overview of the important challenges involving the model validation phase.

Sensing and actuation integration
Energy efficiency, multi-modal operation and complex environments (e.g. unsteady air flow) pose challenges for aerial drones (Floreano & Wood, 2015;Karydis & Kumar, 2017). In this respect a lot about coping with these challenges can be learned from nature (Altshuler et al., 2015;Chin, Matloff, Stowers, Tucci, & Lentink, 2017;Mackenzie, 2012). In nature, distributed mechanoreceptors (e.g. special feathers and hairs) on the wings of birds and bats sense the air flow over the wings (Altshuler et al., 2015;Mohamed et al., 2014). Furthermore, muscles and membranes morph and stiffen the wings of birds and bats Matloff et al., 2020). We expect that sensing and actuation embedding in the material will be essential in being able to i) perceive flow motion and ii) modulate the stiffness of the wings, in order to influence the wing-air interaction, see Fig. 5. This will both be helpful for experimental validation in the experiments described earlier as well as for creating a new robotic bird with unprecedented flight dexterity where stiffness modulation can be used as additional control input.
Extensive research is available on bio-inspired hair flow sensors for all kinds of flow regimes (Jiang et al., 2020;Mark, Xu, & Dickinson, 2019;Mohamed et al., 2014). Specifically, carbon nanotube based hair flow sensors show promising results for the flow regimes of planes (and birds) (Slinker, Kondash, Dickinson, & Baur, 2016). Moreover, the field of flexible endoscopes on the other hand provides extensive research on distributed variable stiffness (Blanc, Delchambre, & Lambert, 2017). Furthermore, distributed variable stiffness for robotics is found in applications with flapping fins, inspired by the fins of fish and dolphins (Huh, Park, & Cho, 2012;Jiang & Gravish, 2018). The mentioned research demonstrates the use of functional materials and structures for flow sensing and variable stiffness mechanisms. For example, variable stiffness fins maximise thrust generation in different modes of locomotion (Jiang & Gravish, 2018) and bio-inspired hair flow sensors predict aerodynamic parameters in turbulent flow (Magar, Pankonien, Reich, & Beblo, 2019). Despite the possibilities that already exist for flow sensing and variable stiffness, both concepts are not properly applied yet in the field of flapping wing aerial robotics. To date, it is for example still unclear how to implement flow sensing arrays in drones due to its complexity (Mark et al., 2019). On the other hand, embedded strain sensors for measuring wing deformation is a more mature technology that can be used in both rigid and flapping wing drones (Shin, Castano, Humbert, & Bergbreiter, 2016;Wissman et al., 2013).
The idea of embedding flow sensing and variable stiffness in the wings fits well with the concept of functional structures. These can be used to change the shape and composition of robots to accommodate opposing dynamic requirements (e.g. both flexible and load-bearing) and to provide additional functionalities (Mintchev & Floreano, 2016). By coupling sensing, actuation and control in one composite material, actual smart or robotic materials can be achieved (McEvoy & Correll, 2015) that use distributed functional materials and structures pushing the state of soft robotics further.
The state-of-the-art in 3D-printing presents promising features for the ambitious fabrication of embedded functionality (McEvoy & Correll, 2015;Rafiee, Farahani, & Therriault, 2020). With the advent of multimaterial 3D-printing it has become possible to generate new multimaterial structures over three or four orders of magnitude in scale, enabling new architectures, embedded electronics or composites directly from CAD models (Hawkes & Cutkosky, 2018;Yang et al., 2018). By means of 3D-printing, electromechanical sensors can directly be embedded into components (Dijkshoorn et al., 2018) and can be used for closed-loop operation (Zolfagharian, Kaynak, & Kouzani, 2020). Furthermore, many of the robots used for studying biological flight are, to a certain extent, manufactured with certain 3D-printing techniques in order to create lightweight, structural components (Ajanic et al., 2020;Chang et al., 2020;Folkertsma et al., 2017;Gerdes et al., 2014). On the other hand the current multi-material 3D-printing methods can only work with a handful of materials, where every material adds complexity to the process and where the various 3D-printing techniques create prints with a limited range of length scales (up to three or four orders of magnitude) (Hawkes & Cutkosky, 2018;Yang et al., 2018). This still limits the achievable resolution and complexity of, for example, embedded sensors (Schouten et al., 2020).

Experimental setups for model validation
Experimental testing plays a crucial role in terms of model validation as well as source of information regarding suitable turbulence models to be used. What is needed in order to validate the dynamic model is a technology able to visualise the flow vector field. Particle Image Velocimetry (PIV) is a class of methods able to measure the flow of a fluid with optical tracking methods, using high speed cameras and intense light sources such as lasers. The resultant velocity vector field can be reconstructed from sufficiently narrow (in time) images of particle clusters. PIV has been used to experimentally verify hypothesised wake patterns (Henningsson, Spedding, & Hedenström, 2008;Spedding, Rosén, & Hedenström, 2003). The use of more complicated PIV setups in order to study unsteady wake patterns connected to flight performance constitutes a stimulating research topic in its own right (Gutierrez, Quinn, Chin, & Lentink, 2017).
The main difficulty in using PIV methods to validate unsteady wake patterns produced by flapping is that a complete 3D fluid flow reconstruction is crucial to validate and understand the model, since flapping flight evolves naturally in 3D. An example is the LEV, a crucial liftenhancing 3D vortex structure on a bird wing that helps a bird land and take-off. This flow structure is commonly attributed to rotational accelerations induced by the flapping motion (Bomphrey, Lawson, Harding, Taylor, & Thomas, 2005;Lentink & Dickinson, 2009). PIV in contrast is typically a 2D method and its extension to the 3D case is highly nontrivial and at the current state comes together with significant drawbacks. Some techniques sacrifice time resolution (Thielicke & Stamhuis, 1994), others significantly limit the measurement volume (Scarano, 2013). Recent advancements in post-processing algorithms for Particle Tracking Velocimetry (PTV) may open the door to large-volume flow velocity analysis (Faleiros, Tuinstra, Sciacchitano, & Scarano, 2019;Schanz, Gesemann, & Schröder, 2016).
Data analysis for PIV methods is a crucial and difficult task for velocity field reconstruction. In addition, at a validation phase, the nature of turbulent flow makes the comparison between measures and model output inherently useless at scales where chaotic patterns emerge. This problem is again a very stimulating one and many contributions are using a wide range of tools to reconstruct the velocity field, from Proper Orthogonal Decomposition (Li, Dong, & Liang, 2016) to deep convolutional neural networks (Lee, Yang, & Yin, 2017).
Another challenge connected to the design of the setup is related to a precise visualisation of the PIV system while measuring the aerodynamic forces on the heaving aerofoil. In fact, in order to study the fluid-solid interaction system with its bilateral dynamic interactions, PIV has to be used to measure flow around an actively controlled mechanical systems. To the knowledge of the authors, current stateof-the-art setups do not allow for these requirements in case of large stroke-length (up to 0.5 m) and flapping frequencies greater than 3 Hz. To overcome these limitations, accurate motion tracking control over the actual flapping profile is required, and this must all be realised in such a manner that the PIV setup still has a clear view of the wing surface throughout the entire stroke cycle. When done correctly, the intensive testing of the heaving aerofoil, will enable the study and isolation of 2D-effects helping and supplementing the study of fully flapping wings and its 3D-effects.
Finally, although wind-tunnel testing is fundamental when it comes to flow-visualisation and data gathering, real-life flight behaviour is not trivially extrapolated from models and wind-tunnel data alone. Eventually, there is no perfect substitute for the actual system flying in open air. Open-air experiments will focus more towards general behaviour such as stability, maneuverability and overall efficiency. A challenging aspect in this context resides in the fact that, due to the unavoidable change in environmental conditions, the comparison between different experiments will not be trivial and will feature statistical analysis and conclusions.
All of the data gathered during wind-tunnel and open-air testing will be used to validate the port-Hamiltonian representation of the fluid and gain information about the turbulence model to be used. Furthermore, with reference to the discussion in Section 5, it will be crucial to tune the control strategies and give rational hints in the sensing and actuation designs.

Conclusions
In this paper we embody the vision of the PortWings project, aiming at pushing the knowledge behind the secrets of flapping flight forward. The core methodology relies on the use of a new modelling paradigm starting from first physical principles, able to deal with the complicated nature of the underlying dynamical system without introducing quasi-steady aerodynamic assumptions. In particular, we extend port-Hamiltonian system theory to a general geometric description of fluid-solid interaction systems, by bringing the Navier-Stokes equation in the pH framework and allowing for a composable description with variable boundaries. The goal is to create a unified theory able to describe flapping dynamics in different regimes, for different maneuvers, reproducing the effects which have been experimentally observed and measured. It is then possible to take advantage of the pH model by developing suitable discretisation methods preserving the power flows of the underlying continuous systems. Reasonably this procedure is able to generate big amounts of physically reliable data corresponding to bird flight regimes. A combination of the use of control paradigms inherited by the pH methodology and novel machine learning methods able to exploit the generated data for optimising actuation strategies, will provide physically based insights (rather than trial and error iterations) for designing the control system able to achieve flight. A discussion on the various technologies that can serve to this purpose is reviewed and discussed. Finally, a robotic bird able to reproduce the different effects and motions used by birds for flying will be realised, towards delivering an undisputable, unique proof of decoding the secrets of flight.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding
This work was supported by the PortWings project funded by the European Research Council [Grant Agreement No. 787675].