Reliability-based design optimization of control parameters for the optimal lift-to-power ratio of clapping-wing micro air vehicles

Abstract

A reliability-based design optimization framework is introduced and applied to optimize control parameters for the optimal lift-to-power ratio of clapping-wing micro air vehicles. Improving this ratio is essential for achieving high load capacity and long endurance in flapping-wing micro air vehicles, especially for clapping-wing types. First, a lift-to-power ratio solver based on smoothed particle hydrodynamics coupled with the finite element method is proposed. To account for the flexible deformation of the wings during high-frequency flapping, an air-solid interaction numerical model is employed to handle quasi-flexible aerodynamics. An experimental platform using particle image velocimetry is established to validate the numerical model and measure power consumption. Then, considering the structural and environmental uncertainties of clapping-wing micro air vehicles measured through experiments, a reliability-based design optimization with an accelerated Kriging model is proposed and solved by particle swarm optimization. Finally, several design cases are examined. Through aerodynamic analysis and reliability-based optimization, optimal control parameters are identified to maximize the lift-to-power ratio of the designed micro air vehicle. The novel method demonstrates rapid convergence, decreasing from 6 to 4 iterations for two control variables and from 9 to 6 iterations for three control variables, thus enabling highly efficient optimization.

Keywords

clapping-wing micro air vehicles reliability-based design optimization air-solid interaction experimental platform accelerated Kriging model particle swarm optimization

Introduction

Micro air vehicles (MAVs) have long been a research focus due to their potential roles in aerial surveillance, target identification, and communication relay.^1,2 Recently, researchers have shown increasing interest in flapping-wing micro air vehicles (FWMAVs). The clap-and-fling mechanism, initially discovered in insect flight, has been widely adopted in FWMAVs’ designs to enhance the lift generation.³ The efficient lift-generation mechanism can help existing MAVs overcome many bottleneck problems in task performance.^4,5 Further research indicates that the double clap-and-fling mechanism is more effective. Using the double clap-and-fling mechanism, clapping-wing micro air vehicles (CWMAVs) have been developed with an X-configuration flapping wing.^5,6 However, CWMAVs still face the problem of excessive power consumption, which is also a challenge in helicopters and multi-rotor unmanned aerial vehicles (UAVs).^7,8 Therefore, it is necessary to optimize the design by adjusting control parameters to maximize the lift-to-power ratio, enabling CWMAVs to play a greater role in various applications.

Some research has been done to optimize the force-to-power ratio. Ünal et al.⁹ found that increasing the lift-to-drag ratio contributes to range and endurance by enhancing fuel efficiency. Li et al.¹⁰ conducted wing parameter optimization experiments to improve lift and efficiency by modifying the layout of the wing vein, aspect ratio, surface area, and flexibility of the leading-edge rod of the CWMAVs. Nguyen et al.¹¹ investigated the effects of wing flexibility on vertical thrust generation and power consumption under hovering conditions for the hovering FWMAV FlowerFly. The thrust-to-power ratios were compared to achieve a wing configuration with high vertical thrust production and low power consumption. Deng et al.¹² explored the aerodynamic performance of CWMAVs using force and flow field measurements to determine the best flapping frequency, aspect ratio, and flexibility. Groen¹³ experimentally evaluated the lift-to-power ratio $(L / P)$ , emphasizing power consumption as a key factor in defining flapping propulsion performance.

However, current research still has some limitations, such as the inability to accurately calculate the aerodynamic characteristics of three-dimensional flexible wings and the lack of consideration of the structural and environmental uncertainties faced by CWMAVs. These simplifications can lead to significant errors in defining the control strategy corresponding to the optimal lift-to-power ratio. This is the main motivation for the present study.

To address the challenge of calculating the lift-to-power ratio, a lift-to-power ratio solver based on smoothed particle hydrodynamics (SPH) combined with the finite element method (FEM) is proposed. Geometric nonlinearities, wing-wing interactions, and unsteady effects in high-frequency flapping motions of CWMAVs make numerical simulation difficult to meet accuracy requirements. The meshless Lagrangian smoothed particle hydrodynamics method, proposed by Monaghan,¹⁴ is well suited to solving these problems, as large fluid deformations do not cause mesh distortion or connectivity issues. Additionally, the SPH method can be coupled with other Lagrangian techniques, such as FEM.^15,16 Such coupling possibilities with finite element methods make SPH a promising solution for aerodynamic modeling of CWMAVs.

To optimize CMWAV control parameters under uncertain conditions, a reliability-based design optimization (RBDO) with an accelerated Kriging model is proposed. Size constraints make it difficult to use industry-standard parts for important components, such as flapping wings, retarding mechanisms, and flapping-wing driving mechanisms. Additionally, metal materials are unsuitable due to weight restrictions, as they have a lower stiffness-to-weight ratio compared to composite materials. Therefore, 3D printing of polyester materials is considered the primary manufacturing approach for CWMAVs, introducing a certain number of uncertainties.¹⁷ Because of these unavoidable uncertainties, the actual aerodynamic performance of CWMAVs always varies around a deterministic value.¹⁸

The flapping frequency $(F)$ , angle of attack (AoA), and forward velocity $(U)$ are crucial and controllable factors that determine the aerodynamic performance of CWMAVs.¹⁹ The uncertainties caused by the mentioned mechanisms and environmental conditions also directly impact these factors, ultimately influencing aerodynamic performance. To evaluate the uncertainties, the concept of reliability is introduced in this study. Reliability is defined as the probability of achieving desired or acceptable performance under uncertainty.²⁰ Traditional reliability analysis methods include the first-order reliability method,²¹ second-order reliability method,²² moment-based method,²³ response surface method,²⁴ and Kriging model.²⁵ The Kriging model offers many advantages due to its high computational efficiency relative to other methods when handling only a few design variables.²⁶ In addition, artificial intelligence-based optimization methods, such as genetic algorithms,²⁷ artificial bee colony algorithms,²⁸ immune plasma algorithms,^29,30 and particle swarm optimization (PSO) algorithms,³¹ have been recognized as advanced methods for solving optimization problems.

This study develops a framework to identify the optimal control parameters for achieving the highest lift-to-power ratio and enhancing flight performance. It validates the effectiveness of the proposed method for both two and three design variables. The novelty of the present study is threefold: (1) a numerical aerodynamic model is proposed for flexible three-dimensional clapping wings with high-frequency motion, which can calculate the lift-to-power ratio to provide accurate input for optimization; (2) an accelerated Kriging model with a moving search window is proposed to build the performance function; (3) a PSO-based RBDO model is proposed to find the optimal flapping frequency, angle of attack, and forward velocity to obtain the highest lift-to-power ratio under uncertainties.

The remainder of this paper is organized as follows. “Methods” section introduces the proposed SPH-FEM model to solve air-solid interaction problems and describes the proposed RBDO to obtain the optimal control strategy. “Results and discussions” section presents and discusses the results. Finally, “Conclusions” section summarizes the conclusions and limitations. “Future work” section presents section recommendations for future research.

Methods

The workflow of the proposed method is detailed in Figure 1, which is divided into two parts: (I) An air-solid interaction model is used to evaluate the aerodynamic performance of the CWMAVs and gather accurate force data. An experiment explores model validation, power consumption calculation, and uncertainty evaluation. (II) A reliability-based lift-to-power optimization model using the accelerated Kriging and PSO optimization methods is established to find the optimal control parameters with high reliability.

Figure 1.

Workflow of the overall procedure.

Three-dimensional flexible clapping wing modeling

The applied CWMAV prototype features a biplane wing design reinforced with carbon rods acting as veins. According to the flight time sequence, the entire flapping duration is divided into six phases, as shown in Figure 2, which was filmed by a high-speed camera. The wingspan is 280 mm, and the weight is 17 g. The flapping mechanism is crucial for achieving the desired flapping motion of the wings.³² Drawing inspiration from the flapping motions of flying creatures, the flapping amplitude is set to 75°. The instroke and outstroke phases indicate the folding and moving apart of the clapping wings, respectively. A clamping regulator is located at the bottom of the CWMAV, which can regulate AoA and provide input for control signals.

Figure 2.

Visualization of a flapping period represented by high-speed camera recordings.

The proposed aerodynamic analysis method (SPH-FEM), which combines SPH and FEM techniques, offers the advantages of being mesh-independent and highly effective at capturing flexible deformation in fluid-structure interaction problems. The three-dimensional flexible clapping wings were innovatively designed by integrating traditional wing geometry parameters to simulate wing deformation during a flapping cycle. A quasi-flexible wing modeling method, as in Figure 3, was adopted to restore the clapping wings, consisting of a combination of ball joints and linear springs that replace the veins of the actual clapping wings. The flexibility was controlled by changing the parameters of the ball joints and the linear springs. Furthermore, to make the material property similar to that of a polyethylene film of the prototype, Young’s modulus of wings was set at 0.2 GPa, Poisson’s ratio at 0.3, the elastic coefficient of restitution at 0.001, and the Coulomb friction coefficient at 0.45. The proposed quasi-flexible modeling method was shown to produce smooth and adaptable wing deformation. The comparison of wing deformation between experiment and simulation is presented in Figure 4, demonstrating a high similarity.

Figure 3.

Quasi-flexible wing modeling method: (a) physical structure of flexible clapping wing and (b) numerical structure of flexible clapping wing.

Figure 4.

Wing deformation between experiment and simulation: (a) clapping wings in experiment and (b) clapping wings in simulation.

SPH-FEM aerodynamic analysis

The open-source code DualSPHysics^33,34 is used to assist in the analysis of CWMAVs aerodynamic performance. The discrete approximation of a physical quantity $β$ for the number of particles $a$ is calculated by

β_{a} = \sum_{b \in Ω} β_{b} W_{a, b} V_{b}

(1)

where $b \in Ω$ denotes the set of neighboring particles and $Ω$ denotes the set of neighboring particles. The kernel function $W_{a, b} = W (| x_{a, b} |, h)$ is calculated as a function of the distance between the particles $(| x_{a, b} | = | x_{a} - x_{b} |)$ and the smoothing length $h$ . The volume of a neighboring particle $j$ is given by $V_{b} = m_{b} / ρ_{b}$ , where $m_{b}$ and $ρ_{b}$ are the mass and density of particle $b$ , respectively.

Mass conservation (continuity) and momentum conservation are defined by

\frac{D ρ}{Dt} + ρ \nabla \cdot u = 0

(2)

\frac{D u}{Dt} = - \frac{1}{ρ} \nabla p + g + Γ

(3)

where $ρ$ is the fluid density, $u = (u, v, w)$ is the velocity vector in the $(x, y, z)$ -directions, respectively, $p$ is the fluid pressure, $g$ is the gravitational acceleration, $Γ$ is the dissipative terms.

A weakly compressible SPH form is adopted as follows:

\frac{d ρ_{a}}{dt} = \sum_{j \in Ω} m_{b} v_{a, b} \cdot \nabla W_{a, b} + δ h c_{0} \sum_{b \in Ω} V_{b} Ψ_{a, b} \cdot \nabla W_{a, b}

(4)

where $v_{a, b} = v_{a} - v_{b}$ represents velocity of the particles. $\nabla W_{a, b}$ are the kernel gradients. $c_{0}$ is the speed of sound. The $δ$ -SPH coefficient, $δ$ , is assumed to be constant. $Ψ_{i, j}$ is the diffusion term.

The air-solid interaction solving framework in this study is based on the multi-body solver “Project Chrono.” DualSPHysics 5.0 interfaces with the rigid-body dynamics solver within Project Chrono’s library, including constraints and collision detection components. Project Chrono creates structure-structure interaction using different contact tracing capabilities to achieve a quasi-flexible wing model.

The implementation of air-solid interaction is shown in Figure 5. The active and reactive forces exerted on bodies are calculated using the position P , angular velocity $Θ$ , and linear velocity V determined by Chrono to update the data for the particles associated with each object.

Figure 5.

Diagram of the Project Chrono air-solid interaction.

Three particle spacings $d_{p}$ of the proposed numerical model $d_{p}$ are compared: 0.008, 0.004, and 0.002 m. Table 1 shows the sensitivity analysis between the experiment and simulation values for different $d_{p}$ values by average lift. It is evident that situations with $d_{p} = 0.004 m$ and $d_{p} = 0.002 m$ have a closer magnitude to the experimental lift. Additionally, comparative analysis shows that setting $d_{p} = 0.002 m$ results in an excessively time-consuming model. When considering the total number of particles, setting $d_{p} = 0.002 m$ results in an excessively time-consuming model, which causes the simulation to take 20 h for every 1 s of simulated flapping time at a compute capability of 8.6 (with NVIDIA GeForce RTX 3060 Ti), which is 10 times longer than setting $d_{p} = 0.004 m$ . Therefore, to strike a balance between sufficient precision and computational expense, a magnitude of $d_{p} = 0.004 m$ is adopted for the parametric study in this study.

Table 1.

Average lift differences between experimental and simulation measurements for different $d_{p}$ (Case 1: $F = 12 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$ ; Case 2: $F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$ ).

No.	Experiment	$d_{p}$ 0.008	$d_{p}$ 0.004	$d_{p}$ 0.002
1	0.1616 N	0.0987 N	0.1601 N	0.1592 N
2	0.1639 N	0.1154 N	0.1666 N	0.1614 N

Table 2 lists the parameters and values adopted for the numerical model.²⁹ The air-solid interaction solution connects wing partitioning and air via a message-passing interface, generating a mesh model with approximately 2,100,000 particles. Figure 3(b) presents the simulation structure of the quasi-flexible clapping wing.

Table 2.

Parameters and values used for the numerical model.

Parameter	Value
SPH kernel	Quintic Wendland
Fluid and solid interaction	Project Chrono
Particle spacing	0.004 m
Smoothing length	$1.2 \sqrt{2} d_{p}$
Density diffusion	0.1
Diffusion term	Fourtakas et al.³⁴
Momentum equation	Laminar viscosity
Particle shifting	Off
Speed of sound	$20 \sqrt{g h_{d}}$
Reference density	$1.34 kg / m^{- 3}$
Polytropic index	7
Viscosity	0.0000148
Time integration	Predictor-corrector
Courant number	0.2

The proposed numerical model of clapping wings focuses on the effects of flapping frequency, angle of attack, and forward velocity on aerodynamic performance. The parameters are illustrated in Figure 6. The maximum force was usually generated at the expense of the maximum power consumption. So, the lift-to-power ratio $(L / P)$ is introduced to evaluate aerodynamic performance instead of lift or power consumption alone. A higher ratio implies more lift generation per unit of power consumption, indicating the superior performance of the clapping-wing motion. Based on the power consumption measured by experiment and the lift calculated by simulation, $L / P$ can be used to analyze the aerodynamic characteristics of CWMAVs.

Figure 6.

Schematic diagram of control variables.

Experimental measurement

The experiment platform shown in Figure 7 aims to provide validation data for the SPH-FEM model, experimental measurements of power consumption, and uncertainty quantification.³⁵ The experiment consists of five parts: measurement, CWMAV prototype, data, environment, and synchronizer. A prototype with a crank rocker mechanism and clapping wings was designed with 17 g. The lift data was obtained by an FC Nano17 sensor (ATI Industrial Automation, $F_{x}, F_{y} = 8 N$ , $F_{z} = 14 N$ , $M_{x}, M_{y}, M_{z} = 0.005 Nm$ , 10.1 g, with a sampling rate of $10000 Hz$ , a resolution of $1 / 160 N$ , and measurement uncertainty of 1%) to determine the accuracy of the proposed numerical models. The comparison of instantaneous lift over a stable flapping period under the condition ( $F = 16 Hz$ , $AoA = 15^{°}$ , and $U = 3.5 m / s$ ) is shown in Figure 8. The curve shows that the simulation and experimental results closely match for the instantaneous lift.

Figure 7.

Diagram of the experiment platform composition.

Figure 8.

Instantaneous lift difference between experimental and simulation measurements ( $F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$ ).

The average lift discrepancies during a whole clapping period between the experimental and SPH-FEM measurements are shown in Table 3. There is a deviation from the experimental result, but the difference in all cases is less than 5%, which is definitely within the acceptable range. The high similarity between the instantaneous lift and the low difference in the average lift indicates that the proposed numerical model meets the accuracy requirements of this study.

Table 3.

Average lift differences between experimental and simulation measurements.

Cases	Experiment (N)	Simulation (N)	Differences (%)
$F = 8 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$	0.1327	0.1356	2.1
$F = 12 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$	0.1616	0.1601	0.9
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 15^{°}$	0.1639	0.1666	1.6
$F = 16 Hz$ , $U = 0.5 m / s$ , $AoA = 15^{°}$	0.0372	0.0385	3.4
$F = 16 Hz$ , $U = 1.5 m / s$ , $AoA = 15^{°}$	0.1080	0.1058	2.0
$F = 16 Hz$ , $U = 2.5 m / s$ , $AoA = 15^{°}$	0.0920	0.0935	1.6
$F = 16 Hz$ , $U = 4.5 m / s$ , $AoA = 15^{°}$	0.2285	0.2302	0.7
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = - 10^{°}$	0.0223	0.0225	0.8
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 0^{°}$	0.0418	0.0413	1.1
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 7 . 5^{°}$	0.0822	0.0841	2.3
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 30^{°}$	0.3047	0.3041	0.2
$F = 16 Hz$ , $U = 3.5 m / s$ , $AoA = 45^{°}$	0.3317	0.3337	0.6

An oscilloscope was assembled to collect current and voltage to calculate the power consumption of the CWMAV, considering the uncertainty of frequency derived from the calculation of the force data period provided by FC-Nano17. The uncertainty of the angle of attack was obtained by calibrating the moving image captured by high-speed cameras. The uncertainty of the forward velocity is obtained directly from the velocity sensors of the wind tunnel. For example, the experimentally obtained probability spectral density (PSD) functions under the flight conditions $F = 12 Hz$ , $AoA = 15^{°}$ , and $U = 3.5 m / s$ are plotted in Figure 9. Under wind tunnel conditions with particle image velocimetry (PIV), the $U$ means the wind velocity. The standard deviations obtained from the PSD function curves are 0.153, 0.462, and 0.05, respectively.

Figure 9.

Probability density distribution for uncertainty quantization of design variables of flapping frequency, angle of attack, and forward velocity.

RBDO

RBDO framework

The aerodynamic performance of the CWMAV is susceptible to random processes arising from unavoidable uncertainties. RBDO in this study aims to solve the optimal control parameters for the best aerodynamic performance while using the experimentally obtained uncertainties. The probability of meeting the desired requirements for CWMAV optimization is as follows:

R = \Pr (G (X) > 0)

(5)

where $R$ is the reliability, $\Pr$ represents the probability function, $G (X)$ is the performance function, and $X$ is a vector containing the design variables that are susceptible to uncertainties.

Compared with traditional design optimization methods, RBDO represents a substantial improvement considering optimization objectives and reliability constraints, as follows:

\begin{matrix} \min cost (d) \\ s . t . {\begin{matrix} R_{j} = \Pr (G_{j} (d) > 0) \geq R_{0} \\ d^{L} \leq d \leq d^{U} \end{matrix} \end{matrix}

(6)

where $cost (d)$ is the cost function of aerodynamic performance, $d$ represents the design variables ( $F$ , $AoA$ , and $U$ ), $R_{j}$ is the $j^{th}$ reliability, $R_{0}$ is the required reliability, $d^{L}$ and $d^{U}$ are the lower and upper boundaries of the design variables, respectively.

Based on the uncertainties measured and analyzed by experiment, the objective of optimization is set as maximizing the lift-to-power ratio, using the $F$ , $AoA$ , and $U$ as design variables. In all cases, the range of design variables is introduced based on the actual flight conditions,¹⁹ as follows:

{\begin{matrix} 4 \leq F \leq 30 \\ - 10 \leq AoA \leq 60 \\ 0 \leq U \leq 6 \end{matrix}

(7)

RBDO for clapping wings using Kriging model

In this study, the aerodynamic characteristics were calculated using the lift obtained from the SPH-FEM simulation and the power consumption obtained from the experimental measurements. Although the proposed SPH-FEM simulations can significantly increase the accuracy and efficiency of lift solving, the computational cost also increases, particularly when multiple design variables are involved. Therefore, the high-fidelity aerodynamic analysis model was approximated using a surrogate model.

For the current optimization problem involving three design variables, Kriging model has proven to be highly effective and accurate.¹⁸ Three Kriging models for different parameters were constructed: the Kriging model for power consumption based on the experimental data, the lift Kriging model based on high-fidelity SPH-FEM simulations, and the lift-to-power ratio Kriging model to combine the data from the first two models. The training strategy includes five steps: 1) Latin hypercube sampling (LHS), 2) analysis using the three initial Kriging models, 3) addition of new samples, 4) error estimation, and 5) search for optimal results using PSO-based RBDO.

Two types of cases with different numbers of variables were discussed in this study. The first case designs two control variables: $d [2] = [F, AoA]$ . The forward velocity $U$ is defined as a constant, considering that the forward velocity is sometimes uncontrollable in CWMAVs control. The second case defines forward velocity as a variable, leading to three variables: $d [3] = [F, AoA, U]$ . The number of initial samples and variables are $N$ and $n$ , respectively, then the initial sample matrix can be expressed as

d^{ini} = [\begin{matrix} d_{1}^{1} & \dots & d_{n}^{1} \\ d_{1}^{2} & ⋱ & d_{n}^{2} \\ ⋮ & ⋱ & ⋮ \\ d_{1}^{N} & \dots & d_{n}^{N} \end{matrix}]

(8)

where the elements in each row represent one round of sampling $d^{i}$ , $L / P^{i}$ is calculated by the lift-to-power ratio Kriging model. The matrix can be rewritten as

D^{ini} = [\begin{matrix} d_{1}^{1} & \dots & d_{n}^{1} & L / P^{1} \\ d_{1}^{2} & ⋱ & d_{n}^{2} & L / P^{2} \\ ⋮ & ⋱ & ⋱ & ⋮ \\ d_{1}^{N} & \dots & d_{n}^{N} & L / P^{N} \end{matrix}]

(9)

New samples were added to train the initial model $D^{ini}$ according to the mean-square error (MSE) between the expected and actual answers, calculated as follows:

MSE = \frac{1}{K} \sum_{k = 1}^{K} [{(\tilde{L / P} (d_{k}) - L / P (d_{k}))}^{2}]

(10)

where $\hat{L / P} (d_{k})$ and $L / P (d_{k})$ are the Kriging prediction and the true responses for the $k^{th}$ sample, respectively. An MSE estimator, shown in Appendix 1, is used to identify the point of maximum MSE, where a new training sample is added. In addition to the maximum MSE point, another maximum MSE point inside a moving search window was regarded as the second added point in the optimization loop. This strategy of adding samples, known as accelerated Kriging, avoids superfluous repetitions when the maximum MSE point lies further from the optimal design location.

The Kriging model with two control variables is used as an example to explain the moving search window mechanism. The moving search window for two design variables can be defined as

\begin{matrix} [d_{i + 1}^{L}, d_{i + 1}^{U}] = [d_{i}^{opt} - 0.5 Q, d_{i}^{opt} + 0.5 Q] \\ Q = q \cdot (d^{L} + d^{U}) \end{matrix}

(11)

where $d^{opt}$ is the optimal result, $Q$ represents the range of the moving search window, and $q$ is a scale factor. In this study, three values of $q$ (q = 0.25, 0.15, and 0.1) are compared.

The initial predictions for six samples are presented in Figure 10(a). Figure 10(b) and (c) show the surface and top views of the MSE, respectively. The accelerated Kriging model is proposed using a moving search window to accelerate convergence. Both the global maximum MSE point and the local maximum MSE point inside a moving search window are added into the loop. The moving search window shown in Figure 10(c) as a white rectangle allows an additional training sample to be closer to the region of interest. When the number of design variables increases to three or more, the moving search window would become a surface or hypersurface constraint.

Figure 10.

Initial accelerated Kriging model for the case with two variables: (a) response, (b) MSE model, and (c) moving search window for choosing new samples.

After adding new training samples, two convergence criteria are established as follows: CC1: The optimal prediction error for two consecutive cycles should be less than a given threshold $c_{1}$ :

\begin{matrix} {Error}^{i} = \frac{| {\tilde{L / P}}_{i}^{opt} - L / P_{i}^{opt} |}{L / P_{i}^{opt}} \\ {\begin{matrix} {Error}^{i} \leq c_{1} \\ {Error}^{i - 1} \leq c_{1} \end{matrix} \end{matrix}

(12)

where $Erro r^{i}$ is the results error of the $i^{th}$ iteration, ${\tilde{L / P}}_{i}^{opt}$ represents prediction at the optimal point, and $L / P_{i}^{opt}$ is the true response.

CC2: The difference between optimal predictions in two consecutive cycles should be less than a given threshold $c_{2}$ :

\frac{| L / P_{i}^{opt} - L / P_{i - 1}^{opt} |}{L / P_{i - 1}^{opt}} \leq c_{2}

(13)

$c_{1} = 2 %$ and $c_{2} = 3 %$ for the design and optimization of CWMAV.

The PSO algorithm was used to accelerate the process of selecting a new training sample. PSO is a known metaheuristic technique inspired by the swarm behavior of flocks of birds and shoals of fish, which enables the search of many design space regions. The setup parameters of the

PSO algorithm adopted in this study are listed in Table 4.

Table 4.

Parameters in the PSO design and optimization.

Parameters	Value
Number of variables	2 or 3
Maximum iterations	1000
Number of swarms	1000
Minimum and maximum inertia weight	0.4 and 0.9
Two acceleration factors	2 and 2

Using the Kriging approximation, equation (6) can be rewritten as

{\begin{matrix} G_{j} (d) = \tilde{L / P} (d) - R_{1} \\ d^{L} \leq d \leq d^{U} \end{matrix}

(14)

where $R_{1}$ is the required magnitude of $L / P$ . The goal is to maximize $L / P$ through the solution of the control parameters under uncertainties. A constraint on the objective was set according to the reliability-based analysis for the cases. The case with two variables requires an $L / P \geq 0.12$ with a reliability of $\geq 0.9974$ that corresponds to the three-sigma concept. For the case with three variables, controllable forward velocity leads to larger $L / P$ . So, the required $L / P$ and reliability of the case with three variables are 0.125 and 0.9974, respectively.

To compare the efficiency of moving windows of different sizes, values of $q = 0.25, 0.15$ , and 0.1 were adopted in the accelerated Kriging for comparative analysis. The results for $q = 0.25$ are presented in Tables 5 and 6 for RBDO using two and three design variables, respectively.

Table 5.

Updating process for two variables and one constant.

Iterations	Design variables $[°, Hz]$	$L / P [N / W]$	MSE	Max MSE	Error [%]	Reliability
0	[12.0882, 14.3255]	0.1303	0.0009	0.0033	3.53	0.9931
1	[10.1847, 20.3851]	0.1310	0.0012	0.0029	3.22	0.9955
2	[9.9595, 20.9481]	0.1319	0.0011	0.0024	2.65	0.9989
3	[10.4473, 21.7624]	0.1306	0.0002	0.0023	1.07	0.9984
4	[10.3579, 21.4935]	0.1301	0.0001	0.0023	1.15	0.9986

Table 6.

Updating process for three variables.

Iterations	Design variables $[°, Hz, m / s]$	$L / P [N / W]$	MSE	Max MSE	Error [%]	Reliability
0	[10.4810, 22.0253, 3.6937]	0.1344	1.83e-4	0.0011	3.96	0.9821
1	[9.87342, 18.2278, 3.9494]	0.1384	2.89e-05	0.0016	3.19	0.9955
2	[9.44062, 18.5966, 3.9392]	0.1391	1.83e-05	0.0015	2.33	0.9978
3	[10.0759, 18.2278, 4.1013]	0.1393	2.82e-05	0.0014	2.12	0.9989
4	[10.8861, 19.1772, 4.1772]	0.1406	4.13e-05	0.0014	1.82	0.9998
5	[10.4810, 18.2278, 4.1772]	0.1414	1.44e-05	0.0013	1.23	0.9995
6	[10.4810, 18.2278, 4.1772]	0.1414	1.18e-05	0.0012	1.23	0.9996

Results and discussions

RBDO with two variables

Due to the instability and uncontrollability of the forward velocity, RBDO with two variables searches for the optimal design of the flapping frequency and angle of attack under constant forward velocity $(U = 3 m / s)$ . Six samples were used to generate the initial Kriging model. In Table 6, for $q =$ 0.25, the optimal design variables are [10.3579, 21.4935] with the highest lift-to-power ratio $L / P = 0.1301$ and the reliability 0.9986, which sacrifice the constraint. Reliability requirements are shown not to be met at the beginning of iterations, but the proposed method can help meet the target requirements in only five iterations. When the moving search window size of accelerated Kriging is decreased to $q = 0.15$ or 0.1, the number of iterations is reduced from 4 to 3.

Although higher $F$ or $AoA$ can result in more lift, the optimization results show that the optimal design point of the highest $L / P$ is not where the lift is at its maximum. This is because higher values of $F$ or $AoA$ can increase power consumption, as validated by experiment. As an extremely lightweight MAV, it cannot use batteries with high energy density or large volume. When the CWMAV needs to carry many sensors and other devices to achieve greater autonomy, a control strategy consisting of $F, AoA, U$ for high $L / P$ must be considered in higher-payload scenarios. For the prototype illustrated in Figure 6, an angle of attack of $10 . 3579^{°}$ and a flapping frequency of $21.4935 Hz$ are optimal design values in Table 6, which can inspire the design concept of the next generation of CWMAVs.

RBDO with three variables

For RBDO involving three variables (F, AoA, and U), 20 initial samples were generated by the LHS method to develop the initial Kriging model. In Table 6, it is found that the optimal point is [10.481,18.2278,4.17722] under $q = 0.25$ . The highest $L / P$ prediction is 0.1414 with a reliability of 0.9996.

Compared to the two-variable case, there are differences in the control parameters for flapping frequency and angle of attack, indicating a strong coupling among these three design variables. Controlling only one or two control parameters is definitely not the optimal control strategy. In addition, a larger magnitude of $L / P (0.1414 \geq 0.1301)$ is achieved in RBDO with three variables. The results suggest that the optimal control strategy concluded for the cases with $F$ and $AoA$ should be adjusted when the forward velocity is changing. In addition, it is found that the optimal design is sensitive to the size of the moving search window, with the number of training iterations decreasing from 6 to 5 or 4.

Effectiveness of the proposed Kriging model

In this study, a moving search window was used to accelerate the convergence of the Kriging model. Compared with the classical Kriging model, the accelerated Kriging model can achieve the same goals with fewer interactions. The higher computational efficiency and accuracy are evident in Table 7.

Table 7.

Comparison between classical and accelerated Kriging models.

Parameter	Classical Kriging model		Accelerated Kriging model
Design variables	2	3	2	3
Optimal design points	[10.3421, 21.3922]	[10.4591, 18.212, 4.1772]	[10.3579, 21.4935]	[10.4810, 18.2278, 4.1772]
Optimal prediction	0.13102	0.1412	0.13009	0.1414
Iterations	6	9	4	6

In addition, a narrower window of $q = 0.15$ or 0.1 can accelerate convergence more than a narrower window of $q = 0.25$ . The narrower moving search window accelerates the search process by reducing the number of iterations and time cost, but a narrower window may cause local convergence for nonlinear problems of a lot of design variables. Therefore, narrower search windows are recommended when the number of design variables is small or for a highly convex design space.

Conclusions

This study first addresses the calculation of the lift-to-power ratio. A high-fidelity aerodynamic analysis model was constructed using the SPH-FEM method to handle the three-dimensional flexible wings of CWMAVs. Experimental data were analyzed to validate the numerical model and measure power consumption. There was a deviation from the experimental result compared to the proposed aerodynamic analysis, but the difference was less than 5% for different flapping frequencies, angles of attack, and forward velocities.

Secondly, given the variety of uncertainties in the flapping flight of CWMAVs, reliability-based design optimization was performed using an accelerated Kriging model to identify optimal key control parameters that maximize the lift-to-power ratio. The uncertainties of the control parameters were quantified by PSD functions using the built experiment. The accelerated Kriging model using a moving search window is proved to be more efficient by the cases with two and three variables. Moreover, a narrower moving search window of the accelerated Kriging model measured by the $q$ value can reduce iterations to speed up the optimization loop. A narrower window of $q = 0.15$ or 0.1 can accelerate convergence more than a narrower window of $q = 0.25$ .

Finally, the model accurately determined the control parameters that yield the maximum lift-to-power ratio under different uncertainties. Through iterations of 4 times, an angle of attack of $10 . 3579^{°}$ and a flapping frequency of 21.4935 Hz are optimal design values in the case of two variables with the highest lift-to-power ratio of $L / P = 0.1301$ and the reliability of $99.86 %$ when the forward velocity is defined as constant $U = 3 m / s$ . Through iterations of 6 times, an angle of attack of $18 . 2278^{°}$ , a flapping frequency of $10.481 Hz$ , and a forward velocity of 4.7722 $m / s$ are optimal design values in the case of three variables with the highest lift-to-power ratio of $L / P = 0.1414$ and the reliability of 99.96%. A larger magnitude of $L / P (0.1414 \geq 0.1301)$ is achieved by adjusting three variables.

In this paper, reliability-based design optimization of control parameters for the optimal lift-to-power ratio of clapping-wing micro air vehicles is novelly proposed, considering different flapping frequencies, angles of attack, and forward velocities. This approach addresses the challenge of achieving high lift in CWMAVs with extremely limited volume.

Future work

This study aims to find the best control parameters to make CWMAVs fly more efficiently. However, aerodynamic analysis shows that the geometry of the wings also significantly impacts lift and power consumption. Therefore, future studies will focus on solving structural parameters using the RBDO framework to maximize the lift-to-power ratio of CWMAVs.

Footnotes

Appendix 1

For an unknown response function $f (x)$ , where $x$ is a vector with the design variables, the response function can be approximated using the Kriging metamodel as

(15)

\hat{f} (x) = h {(x)}^{T} v + Δ (x)

where $\hat{f} (x)$ is the prediction from the Kriging model, $[h_{1} (x), h_{2} (x), \dots, h_{p} (x)]$ , are the vectors of regression function, and $[v_{1}, v_{2}, \dots, v_{p}]$ are the unknown coefficient vectors.

In Eq. 15, $h {(x)}^{T} v$ indicates the prediction trend; $Δ (x)$ is a Gaussian process with zero mean and covariance $cov (Δ (x_{i}), Δ (x_{j}))$ , which is determined as

(16)

cov (Δ (x_{i}), Δ (x_{j})) = σ_{Δ}^{2} R (x_{i}, x_{j})

where $σ_{Δ}^{2}$ is the variance of the Gaussian process and $R (x_{i}, x_{j})$ is the correlation function of the Gaussian process.

For the Kriging meta-model with $n_{is}$ initial samples, $[x_{i}, f (x_{i})]$ , where $i = 1, 2, \dots, n_{is}$ ; the coefficient vector $v$ in Eq. 15 is calculated as follows:

(17)

v = {(H^{T} R^{- 1} H)}^{- 1} H^{T} R^{- 1} f

where $R$ is a correlation matrix with elements $R (x_{i}, x_{j})$ , $i, j = 1, 2, \dots, n_{is}$ , $H = [h {(x_{1})}^{T}, h {(x_{2})}^{T}, \dots, h {(x_{n_{t}})}^{T}]$ and $f = {[f (x_{1}), f (x_{2}), \dots, f (x_{n_{is}})]}^{T}$ .

Using the above equations, the predicted response for a new point $x_{new}$ can then be estimated as follows:

(18)

\hat{f} (x_{new}) = h {(x_{new})}^{T} v + r {(x_{new})}^{T} R^{- 1} (f - Hv)

where

(19)

r (x_{new}) = (R (x_{new}, x_{1}), R (x_{new}, x_{2}), \dots, R (x_{new}, x_{n_{x}}))

The prediction MSE at the new point is then calculated as

(20)

\begin{matrix} MSE (x_{new}) = σ_{Δ}^{2} {1 - r (x_{new}) R^{- 1} r (x_{new}) + \\ {[H^{T} R^{- 1} r (x_{new}) - h (x_{new})]}^{T} {(H^{T} R^{- 1} H)}^{- 1} \\ [H^{T} R^{- 1} r (x_{new}) - h (x_{new})]} \end{matrix}

Handling Editor: Chenhui Liang

ORCID iDs

Yanwei Zhang

Shui Yu

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The data that support the findings of this study are available from the Corresponding Author [Yanwei Zhang], upon reasonable request.*

References

Dong

Zhao

. A balance between aerodynamic and olfactory performance during flight in drosophila. Nat Commun 2018; 9(1): 3215.

Han

Xiang

, et al. Aerodynamic performance and analysis of a hovering micro-scale shrouded rotor in confined environment. Adv Mech Eng 2019; 11(4): 1687814018823327.

Wang

Chen

. Leading-edge vortex enhancement of a flexible flapping wing with the clap-and-fling mechanism. Phys Fluids 2025; 37(1): 011917.

Karásek

Muijres

De Wagter

, et al. A tailless aerial robotic flapper reveals that flies use torque coupling in rapid banked turns. Science 2018; 361(6407): 1089–1094.

Zhang

. Optimized actuator design for flapping-wing robots: a multi-objective approach to mimic natural flapping dynamics. Adv Mech Eng 2025; 17(4): 16878132251335551.

Kajak

Karásek

Chu

, et al. A minimal longitudinal dynamic model of a tailless flapping wing robot for control design. Bioinspir Biomim 2019; 14(4): 046008.

Sal

. Effects of the actively morphing root chord and taper on helicopter energy. Aircr Eng Aerosp Technol 2020; 92(2): 264–270.

Yesilbas

Özgür

Ozen

, et al. Simultaneous and stochastic design of piston-prop tuav vertical tail and its autonomous system. Aircr Eng Aerosp Technol 2025; 97(4): 395–404.

Ünal

Öz

Ünal

, et al. Enhancing aerodynamic performance of a two-dimensional airfoil using plasma actuators. Aerosp Sci Technol 2025; 158: 109882.

10.

Shen

, et al. Experimental study on the wing parameter optimization of flapping-wing aircraft based on the clap-and-fling mechanism. Int J Aeronaut Space Sci 2022; 23(2): 265–276.

11.

Nguyen

Chan

Debiasi

. Experimental investigation of wing flexibility on force generation of a hovering flapping wing micro air vehicle with double wing clap-and-fling effects. Int J Micro Air Veh 2017; 9(3): 187–197.

12.

Deng

Wang

Liu

. Experimental study of a bio-inspired flapping wing MAV by means of force and PIV measurements. Aerosp Sci Technol 2019; 94: 105382.

13.

Groen

. PIV and force measurements on the flapping-wing MAV DelFly II. MSc thesis, Delft University of Technology, 2010.

14.

Monaghan

. Smoothed particle hydrodynamics. Rep Prog Phys 2005; 68(8): 1703–1759.

15.

Fragassa

Topalovic

Pavlovic

, et al. Dealing with the effect of air in fluid structure interaction by coupled SPH-FEM methods. Materials 2019; 12(7): 1162.

16.

Ren

Yan

Yang

, et al. Fast SPH simulation for gaseous fluids. Vis Comput 2016; 32(4): 523–534.

17.

Yin

Kaynak

, et al. Uncertainty and disturbance estimator-based control of a flapping-wing aerial vehicle with unknown backlash-like hysteresis. IEEE Trans Ind Electron 2020; 67(6): 4826–4835.

18.

Dey

Mukhopadhyay

Adhikari

. Metamodel based high-fidelity stochastic analysis of composite laminates: a concise review with critical comparative assessment. Compos Struct 2017; 171: 227–250.

19.

Percin

Oudheusden

BWV

Croon

GCHED

, et al. Force generation and wing deformation characteristics of a flapping-wing micro air vehicle ‘DelFly II’ in hovering flight. Bioinspir Biomim 2016; 11(3): 036014

20.

Wang

, et al. Time-dependent mechanism reliability analysis based on envelope function and vine-copula function. Mech Mach Theory 2019; 134: 667–684.

21.

Zhang

, et al. Time-variant reliability analysis via approximation of the first-crossing PDF. Struct Multidiscipl Optim 2020; 62(5): 2653–2667.

22.

Wang

Mourelatos

, et al. Time-dependent reliability of dynamic systems using subset simulation with splitting over a series of correlated time intervals. J Mech Des 2014; 136(6): 061008.

23.

Yuan

. Parameter estimation for bivariate weibull distribution using generalized moment method for reliability evaluation. Qual Reliab Eng Int 2018; 34(4): 631–640.

24.

Xiao

Zuo

Zhou

. A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis. Reliab Eng Syst Saf 2018; 169: 330–338.

25.

Sun

Wang

, et al. A new kriging-based doe strategy and its application to structural reliability analysis. Adv Mech Eng 2018; 10(3): 1687814018767682.

26.

Yang

Meng

Yang

, et al. Adaptive Kriging-assisted enhanced sparrow search with augmented-Lagrangian first-order reliability method for highly efficient structural reliability analysis. Reliab Eng Syst Saf 2026; 267: 111916.

27.

Momani

Abo-Hammou

Mk. Alsmad

. Solution of inverse kinematics problem using genetic algorithms. Appl Math Inf Sci 2016; 10(1): 225–233.

28.

Gorasso

Wang

. Journal bearing optimization using nonsorted genetic algorithm and artificial bee colony algorithm. Adv Mech Eng 2014; 6: 213548.

29.

Kanehira

Mutsuda

Draycott

, et al. Numerical re-creation of multidirectional waves in a circular basin using a particle based method. Ocean Eng 2020; 209: 107446.

30.

Aslan

Oktay

. Path planning of an unmanned combat aerial vehicle with an extended-treatment-approach-based immune plasma algorithm. Aerospace 2023; 10(5): 487.

31.

Yang

Meng

Keshtegar

, et al. ASVR-GPSO: a novel hybrid active support vector regression and global-best partial swarm optimization for structural reliability analysis. Struct Multidiscipl Optim 2025; 68(10): 196.

32.

Cheng

Zhang

, et al. Design of a hover-capable flapping wing micro air vehicle with abdomen-wing coupled control. Chin J Aeronaut 2026; 39: 103807.

33.

Crespo

AJC

Domínguez

Rogers

, et al. Dualsphysics: open-source parallel CFD solver based on smoothed particle hydrodynamics (SPH). Comput Phys Commun 2015; 187: 204–216.

34.

Fourtakas

Dominguez

Vacondio

, et al. Local uniform stencil (LUST) boundary condition for arbitrary 3-d boundaries in parallel smoothed particle hydrodynamics (SPH) models. Comput Fluids 2019; 190: 346–361.

35.

Fan

Bortoni

Hao

, et al. Upstroke wing clapping in bats and bat-inspired robots offers efficient lift generation. J R Soc Interface 2025; 22(223): 20240590.