{'fraction_non_alphanumeric': 0.020467941995443533, 'fraction_numerical': 0.008678986005135067, 'mean_word_length': 4.981721825654337, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the parametrisation of motion kinematics for experimental aerodynamic optimisation

The levels of agility and flight or swimming performance demonstrated by insects, birds, fish, and even some aquatic invertebrates, are often vastly superior to what even the most advanced human-engineered vehicles operating in the same regimes are capable of. Key to this superior locomotion is the animal's manipulation of the generation and shedding of vortices through optimal control of their motion kinematics. Many research efforts related to biological and bio-inspired propulsion focus on understanding the influence of the motion kinematics on the propulsion performance and on optimising the kinematics to improve efficiency or manoeuvrability. One of the first challenges to tackle when conducting a numerical or experimental optimisation of motion kinematics of objects moving through a fluid is the parameterisation of the motion kinematics. In this paper, we present three different approaches to parameterise kinematics, using a set of control points that are connected by a spline interpolation, a finite Fourier series, and a reduced order modal reconstruction based on a proper orthogonal decomposition of a set of random walk trajectories. We compare the results and performance of the different parameterisations for the example of an experimental multi-objective optimisation of the pitching kinematics of a robotic flapping wing device. The optimisation was conducted using a genetic algorithm with the objective to maximise stroke average lift and efficiency. The performance is evaluated with regard to the diversity of the randomly created initial populations, the convergence behaviour of the optimisation, and the final Pareto fronts with their corresponding fitness values.


Introduction
Unsteady locomotion of flying animals such as birds and insects have inspired the development of micro air vehicles, which are less than 15 cm in size and operate in a Reynolds number range from 10 to 10 000 [21,20,22]. Their field of application includes information gathering in confined spaces, areal mapping, or the transport of small goods [7,1]. Flapping wing configurations are a suitable manner of propulsion for micro air vehicles, since they offer better efficiency at low Reynolds number (Re < 100) and an improved manoeuvrability in comparison to fixed and rotating wing configurations [3,11]. The development of complex flapping wing kinematics for various flight scenarios poses a challenge in this field. A similar challenge is the selection of optimal motion kinematics of other bio-inspired propulsion systems such as underwater vehicles mimicking fish propulsion. A standardised optimisation procedure facilitates the identification of motions with desired characteristics. The motion kinematics have to be parameterised by a finite number of variable parameters prior to optimisation to reduce the dimensionality of the solution space and allow for optimal solutions to be found in an affordable way within a reasonable arXiv:2110.07677v1 [physics.flu-dyn] 14 Oct 2021 timeframe. The optimisation solution landscape can dependent strongly on the parameterisation and possible optima might be missed [17]. The ideal parameterisation has as little variable parameters as necessary to reduce the computational complexity of the optimisation problem without restricting the solution space or creating a bias in a certain direction. Yet, any parameterisation will exclude possible solutions due to discretisation in comparison to a continuous solution in an infinite-dimensional solution space and it is vital to carefully select the most suitable approach for each optimisation problem at hand.
Past efforts to optimise bio-inspired motion kinematics initially focussed on simple linear and harmonic motions. Tuncer and Kaya [26] numerically optimised a flapping airfoil for combined maximum thrust and efficiency. They parameterised the sinusoidal plunge and pitching motions with a variable amplitude and phase shift for a fixed frequency. Berman and Wang [4] implemented a quasi-steady model of insects with a hybrid optimisation algorithm. Their kinematic function enables a continuous transition of the pitching angle evolution β with the parameter C between sinusoidal and trapezoidal motions according to with f being the flapping frequency, Φ the phase shift between the stroke and pitching motions, β 0 the offset angle and a scaling factor β m for further modulations. The temporal evolutions of the stroke and elevation angles were defined similarly, yielding a total of twelve variable parameters. Solutions for these twelve parameters were determined that maximised the aerodynamic efficiency while providing enough lift to support the body weight of different insects. The aerodynamic performance can be further increased by more complex, non-harmonic, and asymmetric kinematics, which have to be defined using more parameters or different base functions. The pitching kinematics of insects in nature are often asymmetric which could be an evolutionary adaption to improve the performance [15]. Martin and Gharib [18] were among the first to implement a kinematic motion function with the possibility to represent asymmetric motions. Liu and Aono [15] performed an optimisation of the motion of a pectoral fish fin which was defined as a complex combination of trigonometric expressions for all three spatial angles. The motion kinematics were successfully optimised for a combined objective of minimal side-thrust and maximal efficiency. They noticed an increasing difficulty in finding the global optimum with increasing number of parameters describing the kinematics. Mandre et al [17] investi-gated the influence of the parameterisation with a finite Fourier series on a heaving and pitching hydrofoil. Small modifications of the parameterisation can lead to a significantly altered solution landscape with a bias towards small but pronounced local optima. Gehrke and Mulleners [8] studied pitching motions with increased complexity, which are able to mimic more closely the motions observed in nature. The stroke motion was fixed and the elevation motion was omitted. The pitching angle was defined by a set of control points in the angle-time domain, which were connected with a third-degree spline. The resulting optimal pitching kinematics forming the Pareto front in the stroke average lift versus efficiency diagram vary significantly and highlight the potential of asymmetric and more complex kinematics to push the performance envelope of bio-inspired micro aerial vehicles.
Here, we will present three different approaches to define and parameterise kinematics for optimisation studies. The different approaches will be presented, compared, and evaluated for the example of the experimental optimisation of the pitching kinematics of a flapping wing, building upon the work presented in Gehrke et al [9] and Gehrke and Mulleners [8]. An optimisation is performed for each kinematic function with 6, 12 and 18 parameters, leading to more than 30 000 individual experiments in total. The different approaches to parameterise motion kinematics for optimisation selected here are: control points connected by a spline interpolation, a finite Fourier series, and a linear combination of modes determined by a modal decomposition of kinematics created by a random walk 2 Materials and Methods

Experimental set-up
The kinematics for the different optimisation studies are evaluated with a robotic flapping wing mechanism immersed in an octagonal water tank (see also [9,8]). The wing is a rectangular flat plate with a chord of c = 34 mm and a span of R = 107 mm ( fig. 1). The flapping frequency of the system is f = 0.25 Hz with a peak-to-peak stroke amplitude of φ = 180°. The flapping frequency and stroke amplitude remain constant throughout this study and lead to a Reynolds number of Re = 2φf cR 2 /ν = 4895 at the radius of the second moment of area R 2 , for a kinematic viscosity of ν 20°C = 1.00 × 10 −6 m 2 /s for all experiments [8,23].
These characteristic parameters are selected to match typical flapping wing vehicles and larger insects such as the hawk moth [15,22]. The wing is actuated by two servo motors (Maxon motors, type RE35, 90 W, 100 N mm torque, Switzerland) at the top of a 2-axis shaft system which controls the pitch and stroke angles of the wing ( fig. 1). The wing moves for a total tip-totip amplitude of 0.47 m which results in a 6.91 c tip clearance to the tank for an outer diameter of 0.75 m.
Previous studies have confirmed that no wall effects are present for a tip clearance larger than 5 c [16,12]. The flapping wing apparatus was specifically designed to run unsupervised for large amounts of time and to execute complex motion profiles with high repeatability. Preliminary tests have determined the maximum error between prescribed motion and motor response on the mechanism to be < 0.1°for high acceleration kinematics [8]. The motors are controlled by a motion control board (DMC-4040, Galil Motion Control, USA). A sixaxis force and torque transducer (Nano17, ATI Industrial Automation, USA) is installed at the wing root in fig. 1b to capture the time resolved aerodynamic loads. The sensor is calibrated with SI-12-0.12 to record at a resolution of 3.13 mN for force and 0.0156 N mm for torque measurements. A data acquisition card (National Instruments, USA) is used to record the force and torque transducer signal at a sample frequency of 1000 Hz. The lift and power data was filtered with a zero phase delay low-pass 5th order digital Butterworth filter for the phase-averaged, time-resolved plots. The filter's cut-off frequency was selected to be 12 times the hovering frequency of the system.

Kinematic functions
Three different approaches to parameterise motion kinematics are selected here: a) control points connected by a spline interpolation, b) a finite Fourier series, and c) a linear combination of modes determined by a modal decomposition of kinematics created by a random walk. The different parameterisation were selected by their ability to mimic trapezoidal and sinusoidal pitch angle profiles commonly seen on robotic flapping wing devices as well as more complex pitch angle kinematics observed on nature's fliers ( fig. 1c). We created kinematic functions using the three approaches with three different parameter counts (p = 6, 12, 18). The first (p − 1) parameters are used to modify the pitching angle evolution and the phase shift is controlled with the last parameter ∆t 0 for all kinematic functions. The motions are described by the pitching angle β ( fig. 1). The sinusoidal stroke angle φ remains unchanged throughout Fig. 1 a. Flapping wing mechanism with stroke angle φ and pitch angle β for the hovering motion, b. wing geometry and axis of rotation, c. pitch angle β kinematics of various natural fliers, adopted from [15] and trapezoidal pitch angle kinematics from a robotic flapper [5].
the optimisations and the elevation angle is kept at zero due to its negligible influence on the lift production [15].

Control points connected by splines
As a first intuitive approach, we select a limited number of control points that will be connected using spline interpolation. To start, we create a symmetric pitch profile by distributing (p − 1) control points equidistantly throughout a half-stroke which represents the most efficient use of the available parameters. Freely spaced control points could lead to large gradients for higher parameter counts which can not be executed by our experimental set-up. In a way, the temporal spacing of the control points acts as a build-in low-pass filter for the pitch angle gradients. If the phase of the control points is not fixed, the phase bounds of the individual points depend on each other which makes the selection of the kinematics slightly more complicated. Depending on the specific problem at hand, it is possible to place more control points in certain parts of the motion where larger variations are required to optimise and reduce the parameter count. The beginning and end of the half-stroke are fixed at β = 0°and the control points are solely specified by their pitch angle β n , with 1 ≤ n ≤ (p − 1) as indicated in fig. 2. Next, the control points are connected by a fifth degree spline interpolation to minimise the local accelerationβ. By point-mirroring the first half-stroke around (t/T, β) = (0.5, 0), we obtain a symmetric pitching profile for the back and forth stroke, represented by the grey curve in Fig. 2 Example of the parameterisation of the pitch angle kinematics using control points connected by a spline interpolation with p = 6 parameters. Five parameters describe the pitch angles (β 1 , ..., β 5 ) of the five control points that are equidistantly distributed within a half stroke. A fifth degree spline connects the five control points and the half half stroke profile is point mirrored to create the full symmetric pitch angle profile in grey. The sixth parameter ∆t 0 introduced a phase shift leading to the final pitching profile in colour. fig. 2. Finally, we added a phase-shifted ∆t 0 to allow for advanced and delayed wing rotation with respect to the stroke reversal, leading to the coloured curve in fig. 2. This procedure allows us to created complex nonlinear and non-harmonic motions with high curvatures. The main advantage of the control point parameterisation, is its intuitive and direct local control of the pitching profiles. Local adjustments can be made by moving individual control points and the values of the control points are directly linked to the pitch angle values, which facilitates the formulation of pitch angle constraints due to mechanical and motor limitations. The distribution of the control points along the time axis can be easily varied to adapt the approach to specific restrictions or requirements for a broad range of applications. On the downside, we noticed that the creation of a random set of control point splines is computationally demanding and required 22.8 s on average for motions with 12 parameters on a typical desktop computer. The computational time increases exponentially with the parameter count and takes a couple of minutes for the tests with 18 parameters.

Fourier-series kinematics
In the second approach, we use a finite Fourier sinecosine series to define our pitching kinematics: with a k the cosine coefficients, b k the sine coefficients, and f 0 the flapping frequency. For the current application, we use only the odd coefficients. The first odd sine-coefficients b k are responsible for the base oscillation whereas the odd cosine-coefficients a k are used for smaller modulations. The ability to represent higher curvatures, as they appear in trapezoidal motions, mainly depends on the number of Fourier terms included. A finite Fourier series is not ideal to describe sudden jumps or discontinuities and tends to display an overshoot followed by decaying oscillations near sharp gradients. This behaviour is known as the Gibbs phenomenon and can be reduced by adding a Lanczos-σ-factor [10]. The Lanczos-σ-factor is defined as: with m as the last plus one summing index of a finite Fourier series. The Fourier series based pitching kinematics used in this paper are described by a finite Fourier series with Lanczos-σ correction factor: with (p − 1) parameters, (p − 1)/2 cosines and (p − 1)/2 sine coefficients. The obtained motion kinematics are periodic, continuous in all derivatives and can easily be expressed in a closed form. The half strokes can be asymmetric around the quarter stroke. The phase-shift parameter ∆t 0 is again applied at the end to create advanced or delayed wing rotations with respect to the stroke reversal. Exemplary motions created using 6, 12, and 18 parameters are presented in fig. 3.

Modal reconstruction
In the third approach, we use a low-order eigenmode reconstruction. The eigenmodes are obtained by a modal decomposition of a family of available arbitrary kinematics. This parameterisation approach is particularly elegant if a family of kinematics is known. For example, if a database of kinematics is available from direct observations of insects or fish motions, we can first reduce the dimensionality of the available kinematics using modal decomposition. Measured and new kinematics can then be created by low-order reconstruction of the most dominant kinematic modes. Here, we do not have a database of measured kinematics at our disposal. To demonstrate the concept of eigenmode reconstruction as a parameterisation approach, we have artificially generated a family of randomised kinematics based on a biased random walk algorithm. The random walk motion generator is not an essential part of the approach it merely serves as bypass to obtain a generalised motion database. The random walk algorithm takes randomly sized steps in a discretised half-stroke. The random expectation value distribution for the individual step sizes is defined such that their sum equals zero over one half-stroke. Each step is defined as the expectation value plus a random fluctuation. The mo-tions starts at β = 0°and is set again to zero at the end of the half-stroke. This creates a wider spread of values and higher gradients at the end of the half-stroke in comparison to the start of the half-stroke. This bias is resolved by line-mirroring motions at mid-half-stroke (T /4). An example of a random walk-based motion is presented in the top left panel of fig. 4. A large amount of motions (O(10 5 )) is created within a few minutes with this method.
The family of random walk-based motions are then decomposed using proper orthogonal modal decomposition (POD). All modes start and end with a zero passage and are nearly symmetric ( fig. 4). The symmetry is expected to improve with a larger data basis of random kinematics whereas the computational complexity of the POD is the limiting factor for the number of input motions. The modes exhibit high similarity with Legendre polynomials starting from the second polynomial. The POD ensures the optimality of the projection space for a given set of training motions and the optimal use of parameters. A large variety of complex, asymmetric motions and oscillations can be defined by a linear combination of the first (p − 1) eigenmodes: with ψ n the POD eigenmodes, λ n the corresponding eigenvalues, and a n the parameterisation coefficients.
The first 17 eigenmodes, ψ 1 -ψ 17 , are presented in fig. 4 for the proper orthogonal decomposition of the family of random walk based motions used in this paper. The coefficients a n are considered normalised and have values confined between −1 and −1. The motions are point-mirrored around (T /2, 0) to obtain a full stroke and the phase shift is applied to create the final motion as depicted in the bottom row of fig. 4. Discontinuities in the gradient between the end and beginning of each half-stroke cannot be executed by the motors and are mitigated by a robust local regression smoothing, which is tuned to mainly affect the sections around stroke reversal. Discontinuities in the velocityβ appear for randomly created kinematics in the early populations and diminish when the convergence progresses.

Genetic algorithm optimisation
The pitching kinematics of our robotic flapping wing device have been optimised using the multi-objective genetic algorithm optimisation algorithm (gamultiobj) from the global optimisation toolbox of Matlab ® [6]. A genetic algorithm is a meta-heuristic optimisation method, suited to find the global optimum in a non-linear solution space with a large number of degrees of freedom. It mimics the process of natural selection known from evolution by recombination and mutation. An initial population of 200 individuals, with a uniform distribution of the parameters within their bounds is created to start the optimisation. Each individual corresponds to one periodic flapping motion. The fitness values for two optimisation objectives are directly measured during experiments for each motion. The optimisation objectives are the maximum stroke average lift C L and the maximum stroke average efficiency η. The stroke average efficiency η is defined here as: The overline indicates stroke average quantities. The lift and power coefficients,C L and C P , are defined as The dimensional lift L and power P are normalised with the density ρ, the wing span R, the chord length c and the stroke average velocity U . Each motion is executed over eight cycles. The stroke average quantities are determined based on the last four cycles to exclude any influence of start-up transients. At the selected Reynolds number, the load responses are extremely repeatable and the difference between the individual cycle averaged lift coefficients and the average coefficient over the last four cycles is below 5 % after the first cycle and be- The experimental set-up is subjected to mechanical and safety constraints to protect the equipment. The maximum pitching angle is limited by the range of free movement in the experimental set up such that |β max | ≤ 90°. The minimum angle of β min = 0°is imposed for the parameterisation of the initial half stroke before applying the phase-shift ∆t 0 to avoid ambiguous solutions where advance or delayed rotation is not governed entirely by the phase-shift. The phase shift is bound by [−π/2, π/2]. The minimum and maximum pitch angle bounds are easy to define for the control point approach by directly limiting the values of β of the control points. This does not apply for the Fourier series and modal reconstruction based parameterisations. Here, the extreme angles are the result of a combination of different parameters. We opted here to keep the allowed param-eter ranges as wide as possible and to reject proposed kinematics that are not safe or not feasible for the experimental mechanism before execution.
The most important constraint for our robotic flapper is the pitch angle acceleration. High accelerations put stress on the load cell due to the inertial forces on the wing and have to be limited to avoid damage. The maximally allowed pitch angle accelerations are |β max | ≤ 1800°/s 2 . The majority of randomly created individuals for the control point and modal reconstruction approach exceed the acceleration constraint. This leads to a validity rate of ≤ 10 % in combination with the constraint |β max | ≤ 90°. The parameter bounds are tightened for the Fourier series approach to keep the share of invalid motions within a range of 10 % to 15 %. A simple restriction of the parameter space does not have the same effect for the splines connecting control points and the modal reconstruction. Here, a pre-selection of kinematics is necessary to obtain a full initial population of 200 individuals. The average and peak acceleration values of the randomly created motions increase with the parameter count and so does the share of randomly generated motions that cannot be executed. The pre-selection process takes a couple minutes for the modal reconstruction with 6 parameters and up to 4 h for the control points connected by splines with 12 parameters on a standard desktop computer. The control point approach with 18 parameters has very tightly spaced control points, which lead to high local accelerations above the limit value but they only appear for a very short time such that their risk of damaging the set-up is deemed low. The constraint of a minimum angle β min ≥ 0°in combination with the acceleration and velocity constraints is too strict for the modal reconstruction approach to create an initial population in a reasonable amount of time. By loosening the lower angular constraint β min , and allowing β min < 0, we obtain more kinematics that have a maximum acceleration below the mechanically allowed limit. All parameterisation approaches required a customised implementation of the constraints and adapted selection of the parameter bounds. The implementation of the experimental constraints was most straight forward for the control point based approach due to the more direct local access to the kinematics.

Convergence criterion
To monitor the progress of the optimisations and to decide when to consider the result to be converged, we have implemented a generational distance measure. The generational distance allows for a quantitative compar-ison between the optimisations of the different kinematic functions possible. It is a common measure to determine optimisation progress in genetic algorithms if no optimal Pareto front as reference is known [2]. The generational distance (GD) is defined by Lamont and Veldhuizen [14] as: with S the set of Pareto front members of the current generation, P the Pareto set members of the previous generation, F the objective function, and p the selected distance norm. The generational distance describes the average distance the Pareto front members are shifted between successive generations and represents the optimisation progress. The individual distances are measured between the point of the current front and the closest point from the previous front. All distance measures for the efficiency η and stroke over lift coefficient C L are normalised by the average range of values for the first ten generations (or all generations if less that ten generations have been tested). The convergence speed is defined as the gradient of the generational distance.
An optimisation is considered converged if the convergence speed ≤ 5 × 10 −4 for two consecutive generations and if the generational distance ≤ 5 × 10 −3 . The changes of the Pareto front are assumed to be minor after the convergence criterion is fulfilled.

Results and discussion
Here, we will analyse and compare the performance of the three selected parameterisation approaches for optimisation applications by the example of a multiobjective optimisation of the pitching kinematics of a flapping wing. The pitching kinematics are parametrised by a) control points connected by a spline interpolation, b) a finite Fourier series, and c) a modal reconstruction of POD eigenmodes. Each of the approaches has been tested with 6, 12, and 18 parameters. This lead to more than 30 000 experimental iterations that have been executed over a period of several weeks. We will first compare how well the different approaches cover the motion space, discuss the optimisation process and convergence, and finally compare the Pareto front solutions and the optimised kinematics.

Coverage of the kinematic solution space
Suitable parameterisation approaches for optimisation studies should be able to represent a high variety of qualitative shapes and achieve good initial coverage of the solution space to not a priori exclude potential optimal solutions. To quantify and compare the ability of the different approaches to cover the solution space, we randomly created an initial population consisting of 200 individuals for the different parameterisations with 12 parameters without any acceleration constraint and without applying the phase-shift. The probability density maps for these initial populations of pitching kinematics are presented in the top row of fig. 5 to give a visual impression of the coverage of the solution space. The density maps for p = 18 do not show significant differences and lead to the same conclusions as for p = 12. The splines connecting control points and the modal reconstruction functions cover almost all achievable angles in the first half-stroke ( fig. 5a1,c1). The Fourier-series have a slightly narrower band of admissible motions around stroke reversal with lower diversity ( fig. 5b1). The fixed phase-locations of the control points lead to the occurrence of regions in the density map with higher local probabilities indicating an inhomogeneous coverage and slight bias towards low angles at the control points.
A quantitative comparison is presented in fig. 6, where the brighter bars indicate the values without the application of the acceleration constraint. The parameters we calculated to characterise the coverage of the solutions space and diversity of the parameterised motions kinematics include the maximum pitch rateβ max and pitch accelerationβ max , the half-stroke-averaged standard deviation of the pitching angle σ(β), the halfstroke-averaged normalised standard deviation of the pitch rate σ(β)/σ(β) and the pitch acceleration σ(β)/σ(β). The Fourier series-based kinematics have the lowest maximum pitch rate and pitch accelerations and a lower diversity based as indicated by lower values of the standard deviations of the pitch angle, rate, and acceleration compared to the control point approach and the modal reconstruction. The highest gradients are measured among the modal reconstruction solutions if the acceleration constraint is not applied. The control point kinematics exhibit the highest diversity closely followed by the modal reconstruction method if we do not apply the acceleration constraint.
If we apply the acceleration constraint to omit solutions that cannot or should not be executed by our experimental device, we obtain the probability density plots presented in the bottom row of fig. 5. The quantitative measures after application of the acceleration constraint are depicted in fig. 6 by the darker bars. Here, the acceleration is limited toβ ≤ 1800°/s 2 .  The diversity of the population is characterised by the maximum pitch rateβ max , maximum pitching accelerationβ max , the half-strokeaveraged standard deviation of the pitching angle σ(β), the half-stroke-averaged normalised standard deviation of the pitch rate σ(β)/σ(β) and the pitching acceleration σ(β)/σ(β), and the correlation-coefficients for the pitch rate r(β) and pitching acceleration r(β).
The coverage achieved by the initial populations is now reduced for all approaches especially near the stroke reversal. Motions with larger angles around stroke reversal are omitted due to the acceleration limit. This affects the control point approach the most. The splines connecting the control points tend to develop characteristic spikes when the distance between control points decreases which leads to accelerations above the allowed limit. The initial population for Fourier series kinematics exhibits a small band of angles with very low diversity around stroke reversal and a poor coverage for smaller β. This is caused by the different restriction of the parameter space for the Fourier series. The best coverage is achieved for all kinematic functions at moderate angles between t/T = 0.125 − 0.375. The modal reconstruction has the highest density at lower angles compared to the others ( fig. 5c1). Unfortunately, this increases the chance that a randomly selected individual does not meet the acceleration constraint and needs to be omitted and replaced when building he initial population. To obtain a set of 200 executable individuals for the initial population within a reasonable amount of time, the lower angular constraint β min is loosened. This explains the non-zero probability for β min < 0 in fig. 5c2.
All kinematic functions show a significant loss in diversity following the application of the acceleration constraint. The modal reconstruction now outperforms the control point approach and the Fourier series in all measures. Most of the diversity of the initial control point approach without constraint is created over the tightly spaced control points with a large spectrum of admissible angles. This allows for a large variety of angles, velocities, and accelerations, but it also leads to higher number of omitted kinematics. This could potentially be reduced by distributing the control points differently or by also considering their time coordinates as variable parameters.
Alternative measures for the diversity of the randomly selected individuals in an initial population are obtained by determining the correlation-coefficients for the pitch rate r(β) and the pitch acceleration r(β). They are calculated as the ensemble average correlation-coefficient r for each individual with the rest of the population. A higher diversity within a population now leads to a lower correlation coefficient r. Overall, the diversity decreases upon application of the acceleration constraint as the correlation values increase. The control point approach and the modal reconstruction show again a higher diversity than the population of Fourier series.
To more quantitatively compare the similarity between different randomly generated kinematic, we borrow ideas here from a network topology approach proposed by Zhang and Small [27] and more recently used e.g. by Tandon and Sujith [25] for comparing the similarity between trajectories. The similarity between kinematics can be quantified using the maximum cross-correlation coefficient R i,j for each pair of kinematics β i and β j . The maximum cross-correlation coefficient is defined as: var (β i (0 : T )) var (β j (τ : (T + τ )) (9) with T the period, cov(.) the covariance, and var(.) the variance operators. By taking the maximum over all possible phase shifts τ , the resulting value only evaluates the shape similarity and is not influence by a phase shift, e.g. introduced by ∆t 0 here. Another metric proposed by Zhang and Small [27] is the mean absolute error D i,j between two curves, which is defined here as: We normalised D i,j by the period T and expressed the result in degree such that D i,j can be interpreted as an average pitch angle amplitude difference between two kinematics. Two identical curves have a maximum cross-correlation coefficient R i,j = 1 and a minimal mean absolute error D i,j = 0. Both quantities evaluate different aspects, R i,j provides a measure for the shape similarity between curves and D i,j provides a measure for the amplitude variations between curves. We use both here to evaluate the degree of similarity within a randomly generated population.
After calculating R i,j and D i,j between all curves with the initial population, we construct two similarity networks. Each kinematic motion is represented by a node in the network. Different nodes are connected to each other if the shape similarity measure R i,j is higher than a predefined threshold R th or if the amplitude dissimilarity measure D i,j is lower than a predefined threshold D th . We then determine the node degree for each node, which corresponds to the number of nodes it is connected to. Kinematics with a low node degree show similarity only to a low number of other kinematics in the population. The cumulative distributions of the node degrees (n d ) among the initial population for the three different parameterisation approaches with p = 12 are presented in fig. 7. The light and dark shading correspond respectively to randomly selected populations without and with the application of the acceleration constraint. The more diverse a population of kinematics the more kinematics we expect with a low degree of connectivity and the larger the coloured area in fig. 7. The distributions in the top row have been obtained for the networks based on the shape similarity measure R i,j and indicate how diverse the population is in terms of the shape of the kinematics. The distributions in the bottom row have been obtained for the networks based on the amplitude difference measure D i,j and indicate how diverse the population is in terms of the amplitudes of the kinematics.
The unconstraint populations for the control point and the modal reconstruction have a high degree of diversity as the majority of the kinematics have a low node degree and show similarities to only a few other kinematics in the population. The Fourier series performs slightly lower than the two others. By construction, the Fourier series have an inherent similarity in shape and it is not entirely unexpected to find more nodes with a higher node degree for this approach. Almost all cumulative distributions in fig. 7 shift to higher node degrees when the acceleration constraint is applied, indicating that the constraint decreases the diversity both in terms of shapes and in terms of amplitudes. The only exception is the shape diversity for the modal reconstruction which remains nearly unaffected by the constraint. The spline population is the most diverse for the unconstrained case and the modal population exhibits the highest diversity for the constrained case.
We have also used these metrics to analyse the ability of the different approaches to create the insect-like kinematics and trapezoidal motions presented in fig. 1. All approaches except the modal reconstruction are ca-pable to mimic the presented kinematics with a minimum correlation coefficient R ≥ 0.993 and a maximum mean absolute error D ≤ 1°. The modal reconstruction does equally well for most of the motions except for the trapezoid with steep edges which is approximated with a larger mean absolute error of D = 2.3°.
Based on all visual and quantitative comparison presented here, the control point method yields the best coverage of the solution space and the largest diversity within an initial population if no acceleration constraint is applied. The modal reconstruction method performs better than the control point method when the constraint is applied. The Fourier series populations are least diverse and have the lowest coverage, but still perform sufficiently well for our optimisation application.

Optimisation progress and convergence
Due the specific implementation of the genetic algorithm for our flapping wing optimisation, we can not replace individual kinematics that are omitted based on the acceleration constraint. This is not a general limitation of the approach and can be overcome for future applications. However, in the current study, the total number of the executable individuals in first generations is significantly lower than the nominal 200 individuals due to the implementation ( fig. 8). This is in particular the case for the control point and modal reconstruction methods with 12 or 18 parameters. The Fourier series method is less affected and also the optimisation with a lower number of control points has pop-  ulations with consistently close to 85 % of executable individuals. The number of executable individuals per generation increases for most of the optimisations but never reaches full population size with exception of the control point approach with 6 parameters. The control point approach with 12 parameters and the modal reconstruction approach have particularly low numbers of executable individuals in the early generations, but the situation improves in the course of the optimisation. The lower number of executable individuals in the populations decrease the diversity which hampers the search for the optimum solutions and delays convergence.
The convergence of the different approaches is analysed based on the generational distance introduced in eq. (8). The variation of the generational distance with the number of executed motions is presented in fig. 9. Surprisingly, the parameter count is not the most important factor that influences the convergence behaviour. The largest influence on the convergence is the parameterisation approach. The optimisations using the same parameterisation form groups with similar convergence speed according to the generational distance metric.
The modal reconstructions exhibit the highest relative generational distance and convergence speed at the very beginning. The modal reconstruction with 18 parameters has the highest initial generational distance of 0.346 but also the modal approach with 12 (GD=0.129) and 6 parameters (GD=0.125) have initial values that are one magnitude larger than the values obtained for the other conditions. The higher initial generational distance can be explained by a larger diversity of the initial population. The convergence speed behaves proportionally to the generational distance such that all optimisations converge in comparable time. The smallest relative generational distance and convergence speed is observed for the Fourier series. The motions of the initial Fourier population are less diverse and closer to the final solution due to the different constriction of the parameter bounds. The optimisation progress is slower because the differences between the motions are less pronounced and consequently harder to determine. On the other side, higher diversity for the initial populations leads to higher convergence speed.
The total convergence time differs by 40 % between the slowest (control point approach with p = 12) and fastest (Fourier series approach with p = 6) optimisation and the majority takes between 2486 (Fourier series approach with p = 12) and 2742 (control point approach with p = 12) iterations (table 1). The optimisations are not deterministic and uncertainties are introduced due to random initialisation, recombination and mutation. The differences in convergence time are small when considering that the computational complexity of the optimisations increases exponentially with the parameter count. The Fourier series optimisations tend to reach convergence slightly faster than the modal and control point approaches.

Pareto fronts and optimal kinematics
The resulting Pareto fronts for all 9 optimisations are presented in fig. 10. The colours indicate the number of optimisation parameters and the symbols indicate the parameterisation approach. The different Pareto fronts have the same shape but there is a shift between them and they reach slightly different extreme values. The modal reconstruction approach with p = 6 finds the most efficient solutions (up to η ≈ 1.4) and the Fourier series approach with p = 12 finds the solutions that generate the highest stroke average lift (up to C L ≈ 2.2). Overall, the optimisations with the highest number of parameters are most limited in finding highly efficient  solutions. The higher number of parameters might lead to more complex kinematics that require more power to achieve the same performance in lift. For our current example, p = 12 seems to provide the best balance across the different parameterisation approaches. The additional complexity provided by 12 optimisation parameters compared to the p = 6 optimisations is beneficial in improving the aerodynamic performance without increasing the power requirements and penalising the efficiency. Increasing the parameter count leads to slower convergence and kinematics with higher harmonics variations. To aid the genetic algorithm to fine tune the parameters, we have conducted an additional optimisation where we have taken the results for p = 6 as a start point for a p = 12 optimisation using the Fourier series parameterisation approach. The results did not lead to significant differences with respect to the regular p = 12 optimisation and are not presented here. Overall, the differences between Pareto fronts are of the same order of magnitude as the inherent experimental variations. The influence of the number of parameters seems to have a stronger influence than the parameterisation approach.
In general, the performance of the three parameterisations is comparable and we can not identify a clear preference of any one of them. The optimal kinematics that make up the Pareto front are included in fig. 10. The columns with the same background colour correspond to optimisations with the same number of parameters. The rows correspond to different parameterisation approaches. The colour of the β(t/T ) curves for one half stroke varies from grey to green, red, or black with increasing efficiency. All Pareto kinematics show a continuous transition from high lift motions to more efficient ones, which indicates that the solutions are converged and we obtained a well developed Pareto front. For the same number of parameters, the Pareto kinematics show similar characteristic features for the three parameterisation approaches but also display some subtle differences that cause the shifts in the Pareto front.
A clear common feature is the preference for advanced rotation for most Pareto kinematics. Advanced rotation means that the wing has rotated past its horizontal orientation (β = 0) before the end of the stroke motion. By doing so, it starts the next stroke with a lower and more favourable angle of attack. A positive influence of the advancement of the rotation for lift production has been observed before on insects and their robotics counterparts [19,24,13]. The rotational advancement generates a higher stroke average lift but comes at a cost and there is a continuous decrease of the phase shift towards less advanced rotation with increasing efficiency in all our optimisations. The pitch rate at the end of the stroke is strikingly similar for all cases. The variation in phase shift between the highest lift generating and most efficient kinematics tends to decrease for optimisations with fewer parameters. This suggests that lift increase can be generated at lower cost by a more complex pitching kinematics during the entire stroke instead of relying of the advanced rotation ( fig. 10).
Another common feature is the increase in the maximum pitch angle with increasing efficiency. A larger pitch angle corresponds to a lower angle of attack with leads to lower drag and more efficient kinematics. The values of the maximum pitch angles are approximately the same for all parameterisations and parameter counts except for the Fourier series, which reach slightly lower maximum values of β.
The main differences between the Pareto optimal kinematics of the different cases are the number of local maxima and their timing. The control point and modal reconstruction approaches with p = 6 show two local maxima around t/T = 0.125 and t/T = 0.375 ( fig. 10). These two maxima are less pronounced or even absent in the Fourier series kinematics which look more like a smooth trapezoidal profile. The two local maxima become more pronounced for the control point approach with p = 12 and only the second one around t/T = 0.375 remains present when p is increased to 18 for this approach. The modal reconstructions for p = 18 are also characterised by a single peak that is located around t/T = 0.32. The Fourier series with p = 12 and p = 18 have a distinct bump at the beginning of the stroke, around t/T = 0.1. Overall, the complexity of the motions increases with the parameter count and different shapes and kinematics can lead to similar Pareto fronts depending on the choice of the parameterisation and the parameter count. This is due to the complex relationship between the flow development and the growth of the leading edge vortex for flapping wings motions [5,8].
In previous work, we identified the shear layer velocity u s at the leading edge as a characteristic parameter to scale the force response of generalised flapping wing pitching kinematics [8]. The shear layer velocity act as a proxy for the feeding rate of the leading edge vortex. Higher angles of attack or lower pitching angles lead to a higher shear layer velocity u s which in turn leads to a faster growth of the leading edge vortex. This is the common feature of high lift generating kinematics. The stronger leading edge vortex leads to higher lift, but reaches its maximum size and circulation earlier in the cycle. The maximum in lift and leading edge circulation is typically reached around mid-stroke for high lift kinematics. After the vortex reaches its maximum circulation, it lifts off of the wing which leads to a drop in lift in second half of the stroke.
Efficient motions exhibit smaller angles of attack with pitching angle maxima around β max = 70°( fig. 10). They create a smaller leading edge vortex which grows continuously during the majority of the stroke and promote drag-minimised lift creation. The power and lift coefficients are less dominated by the stroke motion as the angle of attack is smaller causing the shear layer velocity u s to be relatively constant during mid-stroke.

Conclusion
In this study, we presented three different approaches to parameterise motion kinematics, demonstrated their application, and evaluated their performance by the example of an experimental optimisation of the pitching kinematics of a robotic flapping wing in hover. The pitch angle kinematics in our application are described by a) control points that are connected by a fifth-order spline interpolation, b) a finite Fourier series, and c) a linear combination of modes determined by a modal decomposition of kinematics created by a random walk. Each of the parameterisation approaches is implemented with three different parameter counts: 6, 12 and 18, which leads to a total of nine optimisations and more than 30 000 experimental iterations conducted over a period of several weeks. The performance of the different approaches has been analysed by comparing the diversity of the solutions and coverage of the motion space, the optimisation process and convergence, and the Pareto front solutions and the optimised kinematics.
The coverage of the solution space was qualitatively evaluated based on a probability density map of the random initial population of 200 individuals. For a more quantitative comparison of the coverage and the diversity of the solutions, we introduced a number of diversity measures including the maximum pitch angle, rate, and acceleration, their standard deviations, and correlation coefficients for the pitch and acceleration. To further quantify the diversity between the kinematics within the initial populations, we calculated similarity networks based on a mutual shape similarity measure and an amplitude error. The cumulative distribution of the node degree in a population is an intuitive quantitative measure to compare the diversity of the kinematics based on their shape and amplitude. Based on all qualitative and quantitative measures, the control point kinematics exhibit the highest diversity closely followed by the modal reconstruction method if no experimental constraints are in place. The Fourier series-based kinematics have the lowest maximum pitch rate, pitch acceleration, and a lower diversity compared to the control point and the modal reconstruction approach.
In most experimental applications, the theoretical parameter space cannot be fully explored due to mechanical, electrical, or other constraints for example related to the measurement equipment. In our case, the main limitation is imposed by the sensitivity and measurement range of the load cell and the performance envelope of the motor that controls the pitching motion. These limitations lead to a pitch acceleration constraint that significantly reduces the coverage and the solution diversity for all three parameterisation approaches. The practical implementation of an acceleration constraint is easier for the Fourier series approach than for the two other approaches. The modal reconstruction parameterisation performs best in the various diversity measures, followed by the control point and the Fourier series approaches when the acceleration constraint is in place.
The generational distance criterion is implemented to monitor the progress of the optimisation and as a convergence measure to end the optimisation. All optimisations converged in comparable time which did not depend on the parameter count. In the presented optimisation study, it is desirable to have fewer parameters as more parameters increase complexity without a clear aerodynamic benefit. The subtle changes in the kinematics lead to variations in the loads that are of the same order as the inherent experimental fluctuation of the system. This can slow down the convergence of the genetic algorithm especially in later stages when only minor aerodynamic improvements are achieved.
The resulting Pareto fronts for the different optimisations are similar in shape and in the range of values of the two objective functions they cover. The parameter count has a stronger influence on the final results than the parameterisation approach. This is encouraging for the robustness of the three parameterisation approaches presented here. The optimal kinematics corresponding to the Pareto fronts have slightly different shapes depending on the parameter count and the parameterisation approach selected but they yield comparable fitness values. The most prominent common features of the kinematics are their shared preference for advanced rotation and the increase in the maximum pitch angle with increasing efficiency. The main differences between the optimal kinematics of the different cases are the number of local maxima and their timing.
The goal of this paper was to present three different methods to parameterise complex motion kinematics for optimisation studies and to compare their performance on an experimental flapping wing system. Each optimisation application has different objective functions, constraints, and specific challenges. There is no one solution that fits all and the different approaches to parameterise the solutions have their specific advantages and disadvantages outlined in table 2. The Fourier series approach is the easiest method to implement but the introduction of parameter constraints and the parameter bounds cannot be as directly contained as for the control point approach. Nevertheless, the Fourier series has an analytical definition and continuous high-order derivatives which make it a promising parameterisation for structurally sensitive applications like vertical-axis wind turbines or flapping-foil energy harvesters. By definition, the control points are more intuitive to constrain and their parameter values can be interpreted directly in terms of their influence on the temporal evolution and amplitude of the kinematics. With the ability to create sharp and zero gradient curves, the control point method is most suitable for intermittent kinematics like burst and coast swimming of fish-and squid-like devices, and flapping wing flight. The modal reconstruction is extremely versatile and adaptable. The initial family of solutions can be generated by a random walk method as presented here, but can also be a family of kinematics measured in nature or in the lab. Building a parameterisation upon an existing library of kinematics can be an extremely powerful tool. Complex animal-like locomotions can be described by only a few modes which would require many parameters being constructed by analytical kinematic definitions.
Depending on the individual complexity of the mode shapes, it is almost impossible to restrict parameters bounds to meet experimental constraints. Instead, motions have to be omitted after generation if they do not meet constrains such as start and end position and acceleration. This can severely hamper the creation of initial populations for a genetic algorithm optimisation, for example. In the presented application of an experimental optimisation of a robotic flapping wing mechanism, all three approaches yielded similar results in a comparable amount of time. Other applications might have different challenges, but we expect that the main