Nonlinear Dyn

https://doi.org/10.1007/s11071-019-05085-4

ORIGINAL PAPER

On higher-order averaging of time-periodic systems:

reconciliation of two averaging techniques
Marco Maggia · Sameh A. Eisa · Haithem E. Taha

Received: 29 November 2018 / Accepted: 17 June 2019

© Springer Nature B.V. 2019

Abstract In this paper we show how higher-order Time-periodic systems · Vibrational stabilization ·
averaging can be used to remedy serious technical Vibrational control
issues with the direct application of the averaging the-
orem. While doing so, we reconcile two higher-order
averaging methodologies that were developed indepen- 1 Introduction
dently using different tools and within different com-
munities: (i) perturbation theory using a near-identity Nature is teeming with periodically forced creatures
transformation and (ii) chronological calculus using (e.g., insects, birds, fish). With the current active
Lie algebraic tools. We provide the underpinning con- research trend of mimicking nature’s designs, there
cepts behind each averaging approach and provide a are plenty of bio-inspired engineering systems that are
mathematical proof for their equivalence up to the represented by nonlinear time-periodic (NLTP) mod-
fourth order. Moreover, we provide a higher-order aver- els such as swimming robots [1–4] and flapping-wing
aging study and analysis for two applications: the clas- micro-air-vehicles [5,6]. Additionally, high-frequency
sical problem of the Kapitza pendulum and the modern periodic control is a common approach for stabilization
application of flapping flight dynamics of micro-air- and path planning for under-actuated mechanical sys-
vehicles and/or insects. tems (with control inputs less than degrees of freedom)
[2,7–9]. Vibrational control is an open-loop stabiliza-
Keywords Averaging theory · Higher-order aver- tion technique via the application of sufficient high-
aging · Chronological calculus · Lie transform · frequency, high-amplitude, periodic inputs. It was first
introduced by Meerkov [10] for linear time-varying
This work was supported by the National Science Foundation, systems. It is well known for its robustness [11] and
Grant CMMI-1709746 and its continuation CMMI-1846308. elegance (stabilization without feedback). By defini-
M. Maggia · S. A. Eisa · H. E. Taha (B)
tion, vibrationally controlled systems are time-periodic
Department of Mechanical and Aerospace Engineering, systems.
University of California, Irvine, Irvine, CA 92697, USA Analyzing the dynamics of these time-varying sys-
H. E. Taha tems is quite more challenging than that of autonomous
e-mail: hetaha@uci.edu (time-invariant) ones. For example, a NLTP system
M. Maggia usually attains equilibrium when all of its states
e-mail: mmaggia@uci.edu undergo periodic solutions. That is, its equilibria may
S. A. Eisa not be represented by fixed points but rather by periodic
e-mail: seisa@uci.edu orbits. Moreover, even when the system is linear and

123
 M. Maggia et al.

NLTP
cal implementation of the Floquet theorem precludes
Model scrutinizing the dynamical behavior of the system on
an analytical level. The second approach is to use the
averaging theorem to transform the NLTP system into a
nonlinear, time-invariant (NLTI), i.e., an autonomous,
Compute system in which the periodic orbit of the original sys-
Averaging Periodic Orbit tem corresponds to a fixed point. Then, the averaging
then Linearize
theorem guarantees that, for high enough frequency,
exponential stability of this fixed point implies expo-
NLTI LTP nential stability of the corresponding periodic orbit.
Model Model Thus, the averaging approach dodges the determination
of the periodic orbit under study and, hence, provides
Apply a compact, constructive technique to ensure the peri-
Find Fixed Point Floquet odic orbit and to analyze its stability (by ensuring the
then Linearize Decomposition corresponding fixed point and studying its stability). In
fact, averaging has been the standard and most common
technique to design vibrational control systems [7,20].
LTI Despite its convenience to analyze NLTP systems,
Model the application of averaging should be performed care-
fully. In fact, there might not be a more expressive
statement about the subtleties of averaging than that by
Spectral Sanders and Verhulst in their book on the topic [17]:
Analysis
To many physicists and astronomers averaging
seems to be such a natural procedure that they do
Fig. 1 The two main approaches used to analyze NLTP systems
not even bother to justify the process. However,
it is important to have a rigorous approximation
theory, since it is precisely the fact that averaging
seems so natural that obscures the pitfalls and
has a fixed point (ẋ = [A(t)]x), its stability cannot be
restrictions of the method.
assessed through the eigenvalues of the system matrix
[A(t)]. That is, in contrast to linear time-invariant sys- Sticking to the standard form of averaging, its direct
tems ẋ = [A]x where a Hurwitz matrix (its eigenval- application may not be admissible and/or practically
ues have strictly negative real parts) implies exponen- useful in many cases. First, the averaging theorem
tial stability, if all the eigenvalues of the time-varying necessitates a special structure (weakly forced sys-
matrix [A(t)] lie in the open left half of the complex tems). Therefore, it cannot be applied directly to
plane for all times, the system may still be unstable; vibrational control systems, which are typically high-
Markus and Yamabe [12] provided a counterexample. amplitude, periodically forced systems. The standard
In general, the stability analysis of NLTP systems remedy for this case is to apply the nonlinear variation
is conducted using one of two main approaches: a of constants formula [7,21] to transform the multi-scale
numerical approach based on the Floquet theorem [13– system into two companion systems, each of which is
15] and an analytical approach based on the averag- amenable to the averaging theorem (see [7,8,22,23]).
ing theorem [16–19]. As shown in Fig. 1, the first It is interesting to note that the use of the variation of
approach requires determination of the periodic orbit constants (VOC) formula before averaging goes back
(i.e., solving the dynamic equations). Then, the Flo- to Lagrange in his analysis of the perturbed two-body
quet theorem is used to analyze the stability of the lin- problem (see [24, pp. 181–184]). However, the VOC
ear, time-periodic (LTP) system obtained by lineariz- formula, which requires computation of flow maps (i.e.,
ing the dynamics about the concerned periodic orbit. solution of part of the differential equation), may not
For non-trivial applications, this approach cannot be be analytically tractable for complex systems, hence
performed analytically. Hence, the inevitable numeri- frittering away the analytical advantage of the averag-

123
Reconciliation of two higher-order averaging techniques

ing analysis. Second, even when the system structure colleagues Agrachev and Gamkrelidze [28,29] devel-
permits analytical application of the VOC (e.g., sim- oped a new branch of calculus to study time-varying
ple mechanical systems [7]), a major issue with the dynamical systems: the chronological calculus. The
averaging theorem is that it is applicable only for high main objective was to provide an exponential repre-
enough oscillation frequency. That is, it only guaran- sentation of the flow along a time-varying vector field,
tees the existence of a threshold frequency ω∗ above in a way that resembles the exponential representation
which stability of the averaged dynamics implies sta- of the flow of autonomous vector fields given by the
bility of the corresponding NLTP system. The predica- Fliess functional expansion [30]. Utilizing these math-
ment is that almost none of the averaging theorems tell ematical tools, Sarychev [31] and Vela [27] developed
how much ω∗ is. As such, vibrational control is typi- higher-order averaging techniques for time-periodic
cally achieved by ensuring that the averaged dynamics systems, generalizing the classical averaging theorem
is stable and then sweeping the frequency (e.g., by sim- to cases where the excitation frequency is not high
ulation) until hopefully the threshold is hit. Then, the enough. Inspired by the Floquet theorem, the main
averaging theorems guarantee stability for the origi- objective was to determine a time-invariant (averaged)
nal time-periodic system for all frequencies above this system whose flow after one period matches the flow
threshold. The main issue stressed in the recent work of the time-periodic system after one period. Hence,
of Berg and Wickramasinghe [20] is that finding a sta- the stability characteristics of the NLTP system can
ble response of the time-periodic system at a certain be deduced from the stability characteristics of the
frequency does not necessarily mean that the threshold averaged system, thereby extending the Floquet the-
is reached. That is, the time-periodic system may be orem to nonlinear systems. Because this approach is
stable at some ω1 and unstable at ω2 > ω1 (clearly, developed using differential geometric techniques, the
ω1 < ω2 < ω∗ ). This finding also points to the pos- higher-order averaged dynamics are given in terms of
sibility for vibrational stabilization at a frequency less Lie brackets between the time-periodic vector field and
than the threshold value, hence relaxing the averaging its integrals.
unfeasible requirements. Throughout several discussions between the third
In this paper, we propose higher-order averaging author and Prof. Ali Nayfeh, the latter always sug-
techniques as potential remedy for the two issues raised gested that the two averaging approaches are just the
above regarding direct (first-order) averaging. In par- same. Yet, it is not clear in the literature whether such
ticular, we show how higher-order averaging captures approaches are equivalent or not. Therefore, we find it
more dynamical features that are typically neglected a good opportunity to reconcile the two higher-order
by direct averaging. Moreover, the threshold frequency averaging techniques in this special issue in the mem-
may decrease as the order of averaging increases. While ory of Prof. Nayfeh, particularly because he made
doing this, we find it necessary to reconcile differ- several important contributions to the first classical
ent averaging approaches in the literature. Two higher- approach [13–15,32].
order averaging techniques can be found in the litera- The main contribution of this paper is to recon-
ture, which are developed separately employing differ- cile these two averaging techniques, under a unified
ent tools. The first technique is based on the classical notation, providing explicit formulas for higher-order
averaging theory developed by Krilov et al. [25,26]; averaging terms. Moreover, we prove that the two
some authors such as Vela [27] refer to it as the approaches are equivalent up to the fourth order. It is
KBM method. This classical approach has its roots noteworthy to mention that, unlike the classical averag-
in the perturbation theory; the main objective is to ing approach, the one based on chronological calculus
find a near-identity transformation that transforms the does not provide a recursive formula for the compu-
time-periodic system into a time-invariant system after tation of higher-order terms. Hence, the equivalence
neglecting terms of arbitrary small order of magnitude. between the two averaging techniques cannot be sys-
The other higher-order averaging technique is based tematically carried out for an arbitrarily high order. In
on the so-called chronological calculus, exploiting this paper, the equivalence is mathematically proven
Lie algebraic tools. In honor of the 70th birthday of only up to the fourth order. Moreover, we show a
the great Russian mathematician Lev Pontryagin, the “conceptual” equivalence between the two approaches
founder of the optimal control theory, his students and which suggests that the mathematical equivalence can

123
 M. Maggia et al.

be extended for terms beyond the fourth order. This ẋ =  f(x, t, ), (1)
equivalence will allow each approach to benefit from where x ∈ Rn , ẋ denotes the derivative of x with respect
tools and results developed in the other branch. After to t, and  is a small parameter such that 0 <   1.
proving such an equivalence, we apply these tech- The vector field f is smooth in x, analytic in , and T -
niques to two applications: (i) a classical application of periodic in t. The averaged system corresponding to
the Kapitza pendulum (inverted pendulum subjected to (1) is
vertical oscillation) and (ii) a modern application on the
ẋ =  f(x) (2)
hovering flight of insects and flapping-wing micro-air- T
vehicles. In these examples, we show how higher-order where f(x) = 1
T 0 f(x, τ, 0) dτ .
averaging captures more dynamical features than direct 1. If x(0) − x(0) = O(), then there exist b,  ∗ ∈ R+
averaging and how it can mitigate the technical issues such that x(t) − x(t) = O() for all t ∈ [0, b/]
associated with the direct application of the averaging and for all  ∈ (0,  ∗ ).
theorem. 2. If the origin x = 0 is an exponentially stable equi-
The paper is organized as follows: In Sect. 2, the librium of system (2) and if x(0) − x(0) = O(),
two approaches (classical averaging and chronological then there exists an  ∗ such that x(t)−x(t) = O()
calculus) are described separately. The reconciliation for all t ≥ 0 and for all  ∈ (0,  ∗ ). Furthermore,
of the two averaging approaches is given in detail in system (1) has a unique, exponentially stable, T -
Sect. 3. Lastly, Sect. 4 discusses the two applications periodic solution xT (t) such that x T (t) ≤ k for
mentioned above. some k ∈ R+ .
2 Two higher-order averaging techniques Sometimes, the first-order averaging approximation
is too poor and the averaged system (2) fails to cap-
NLTP systems are usually characterized by two ture some important characteristics of the NLTP sys-
timescales: a fast timescale and a slow timescale. The tem. In particular, Theorem 1 is valid only for weakly
goal of averaging is to filter the slow behavior by cap- forced systems or for systems with large separation
turing the mean effect of the dynamics over a period of between the two time scales (e.g., the ratio of the body’s
the fast timescale. During a period, the slowly vary- flight natural frequency to the flapping frequency for
ing states (that do not explicitly depend on the fast the hummingbird,  ≈ ωn /ω,   1). Also, no infor-
time) are approximated as constants. As an example, mation is given about how large this separation should
let us consider a hummingbird flapping its wings. The be for the theorem to be valid. For these reasons, higher-
wings move periodically at a frequency ω that is much order averaging terms are introduced. The higher-order
higher than the natural frequency associated with the terms are also time-invariant and are scaled by powers
bird’s body flight ωn (i.e., ω  ωn ). During one flap- of the small parameter   1. In the averaging the-
ping cycle, the body oscillates in all directions. Nev- ory, an r th-order averaged system implies neglecting
ertheless, on the average, it is hovering in the same terms of order  r . As such, this truncation becomes
place (e.g., over a flower). Moreover, the oscillation of more accurate as the order increases and/or the param-
the body around the mean position/attitude is minimal. eter  decreases.
Therefore, averaging seems so natural to analyze such In this section we briefly summarize the two higher-
a dynamical system. Analyzing the averaged dynamics order averaging approaches that are found in the liter-
can be extremely beneficial because the complexity of ature. We start with the classical averaging approach
the system can be noticeably reduced. In particular, the [15,17,18,33] and proceed with the description of the
explicit dependence on time is eliminated and the aver- less known approach based on chronological calculus
aged system becomes autonomous (i.e., it is fast time- [27,31,34].
invariant). The approximation by averaging is justified
by the following theorem, better known simply as the 2.1 Higher-order averaging via a near-identity
averaging theorem. transformation (lie transform)
Theorem 1 Averaging Theorem (Theorem 10.4 in
[16]). Consider the following NLTP system written in This higher-order averaging approach utilizes a partic-
averaging-canonical form ular change of coordinates to approximate the NLTP

123
Reconciliation of two higher-order averaging techniques

system by a NLTI one. The original NLTP system will 2.1.1 First-order averaging
be at times referred to as the exact system. The method-
ology has been developed from perturbation theory (see Let us rewrite Eqs. (1), (5) and (6) as
[13,15,17,18]).
ẋ = f1 (x, t) +  2 f̂(x, t, ), (7)
We begin our analysis by considering a NLTP
dynamical system written in the averaging-canonical x = y + w1 (y, t) +  ŵ(y, t, ),
2
(8)
form1 (1). Since f is smooth in x and analytical in , it ẏ = g1 (y, t) +  2 ĝ(y, t, ), (9)
can be expanded in power series of  as where f̂, ĝ and ŵ denote remainders. Differentiating
∞
  m−1 Eq. (8) with respect to t and substituting by Eq. (9) in
f(x, t, ) = fm (x, t) = f1 (x, t)
m! the resulting equation, we obtain
m=1
 2 ẋ = g1 (y, t) +  2 ĝ(y, t, )
  ∂w 
+ f2 (x, t) + f3 (x, t) + · · ·
2! 3!
(3) ∂w1  1
+ g1 (y, t) +  ĝ(y, t, ) +
2
where the functions fm are computed as, ∂y ∂t

∂ m−1 f  + O( 2 ). (10)
fm (x, t) = m m−1  . (4)
∂ =0 Substituting Eq. (8) into the argument of f1 in Eq. (7),
The objective is to find a transformation that con- we obtain
verts the x-dynamics into a time-invariant system in the
transformed coordinates y. Given a transformation w, ẋ = f1 (y + w1 (y, t), t) + O( 2 ). (11)
we perform the following change of coordinates x → y Equating the terms in Eqs. (10) and (11) of the same
x = y + w(y, t, ), (5) -power, the  1 equality implies
where y ∈ Rn is the transformed state. For  = 0, ∂w1 (y, t)
g1 (y, t) = f1 (y, t) − . (12)
Eq. (5) reduces to the identity x = y; because  is small, ∂t
the transformation in (5) is a quasi-identity, hence the The vector field g1 (y, t) in Eq. (12) is still time-varying.
name near-identity transformation. This transforma- The only way to make g1 time-independent, while w
tion is sometimes referred to as a Lie transform. When is T -periodic, is to take g1 equal to the average, over a
working with NLTP systems, the function w is also period T , of the vector field f1 [17,18]
required to be T -periodic in t. Furthermore, we assume 1 T
w to be smooth in y and analytic in , so we can expand g1 (y) = f 1 (y) = f1 (y, t) dt. (13)
T 0
it in an asymptotic power series of  as
∞ Combining Eqs. (12) and (13) we find an expression
  m−1
w(y, t, ) = wm (y, t) for the first term of the expansion of the transformation
m! w(y, t) as
m=1
When transformation (5) is applied to system (1), we t
obtain a new nonlinear smooth system in the form w1 (y, t) = (f1 (y, τ ) − g1 (y)) dτ + C1 (y), (14)
0
∞
 m where C1 (y) is an arbitrary t-independent and y-
ẏ = gm (y, t), (6)
m! dependent constant of integration.
m=1
where the functions gm (y, t) depend on the particular
choice of the near-identity transformation (5) and are 2.1.2 Higher-order averaging
in general time-varying vector fields. We now demon-
strate the procedure to make g1 (y, t) time-invariant, The previous reasoning can be easily extended to higher
and we then extend the reasoning to higher-order terms orders. In fact, the terms gm (y, t) are obtained by differ-
gm (y, t) with m > 1. entiating both sides of Eq. (5) with respect to time, sub-
stituting the -expanded Eqs. (1) and (6) in the resulting
1 It is important to notice that, if direct averaging is applied to expression, and finally equating terms of like power of
a system that is not written in the averaging-canonical form, the
following analysis might not be correct. Averaging a system not
. We thus obtain expressions equivalent to Eq. (12)
written in canonical form is known as crude averaging and is (m) ∂wm (y, t)
discouraged as pointed out in [17]. gm (y, t) = f0 (y, t) − . (15)
∂t

123
 M. Maggia et al.

Equation (15) is also known as the homological equa- To compute a certain gm , one only needs to know
tion [17]. Proceeding analogously, we make gm (y, t) the power expansions f1 , f2 , . . . , fm in Eq. (4) and
time-invariant by taking to set the integration constants C1 , C2 , . . . , Cm−1 in
(m)1 T (m) Algorithm 1. The constants of integration Cm (y) are
gm (y) = f 0 (y) = f (y, t) dt, (16) independent of time, but are dependent on the state y;
T 0 0
therefore, different choices of such constants can lead
and the higher-order terms of the transformation w are
to significantly different results. The selection of Cm
given, similarly to Eq. (14) as
is arbitrary; there are typically three main choices that
t 
wm (y, t) = f0(m) (y, τ ) − gm (y) dτ + Cm (y). are found in the literature:
0
1. Stroboscopic Averaging:
(17)
It is important to highlight that in Eqs. (15), (16) Cm = 0. (19)
(m)
and (17), the functions f0 (y, t) are not to be con- By choosing the stroboscopic averaging, the near-
fused with the functions fm (y, t) previously defined identity transformation (5) reduces to the identity
in Eq. (4). In fact, as previously described, the func- when t = nT , with n ∈ Z. In other words, if one
(m)
tions f0 (y, t) stem from equating coefficients of like were to observe, using a stroboscopic light flashing
power of  and can be computed recursively using the at these times, the solutions of the exact system (1)
following algorithm. Given f0(1) = f1(0) = f1 and g1 and the transformed system (21) (with the same ini-
computed as in Eq. (13), then the higher-order aver- tial conditions), they would appear to coincide. For
aged vector fields gm for m > 1 can be computed elegance and simplicity of computations, this is the
following Algorithm 1. In Algorithm 1, we dropped most common choice for higher-order averaging.
the argument of f, g and w for the sake of notation 2. Zero-mean wm :
clarity. Moreover, Aik = k!/((k − i)!i!) are the bino-  τ 
1 T
mial coefficients and [·, ·] represents the Lie bracket Cm = − f̃m dτ dt. (20)
between vector fields. The standard results in the liter- T 0 0
ature [15,33] are written in this section in terms of Lie By choosing the constants of integration as in
brackets to facilitate the reconciliation with the other Eq. (20), the periodic functions wm will have
averaging approach based on chronological calculus. zero
 mean over the period T . In other words,
1 T
The Lie bracket [X, Y] between the two vector fields T 0 wm dt = 0. This choice is preferable when
X(x, t), Y(x, t) ∈ Rn is defined as (see [21,35–37]) fm are written as Fourier series so that wm maintain
∂Y ∂X the same structure (see [18, Ch. 6]).
[X, Y] = X− Y. (18) 3. Reduced Averaging: Cm are chosen so that all the
∂x ∂x
functions gm vanish for m > 1. In other words, the
Algorithm 1 Computing gm for m > 1. dynamics is not captured by the averaged system,
For k = 2, . . . , m: but only via the near-identity transformation. (see
[17, Ch. 3]).
(0)
fk = fk ,
The transformed system (6), when all the vector
t  
(k−1) fields gm are computed as in Eq. (13) and Algorithm 1
wk−1 = f0 − gk−1 dτ + Ck−1 ,
0 (i.e., gm are all time-invariant), becomes a NLTI system

k−2
(1)
fk−1 = fk +
(0) (0)
Aik−1 wi+1 , fk−i−1 , ∞
 m
i=0 ẏ = gm (y). (21)
m!
( j) ( j−1)

k− j
k− j
m=1
fk− j = fk− j+1 + Ai
After proving equivalence with the complete averaging
i=0
( j−1) (k−i−1)
of Sarychev [38] and Vela [27] in the next section, it will
× wi+1 , fk− j−i − f0 + gk−i−1 , j = 2, . . . , k, be clearly seen that system (21) represents the complete
1 T averaging of the NLTP system (1). That is, system (21)
(k)
gk = f0 dt.
T 0 is an exact transformation of the original system (1):
No approximation has been introduced yet; if y(t) is

123
Reconciliation of two higher-order averaging techniques

the solution of system (21), one can exactly recover where

the solution x(t) of the original system (1) via the Lie t
transformation (5). However, when the summation is w1 (y, t) = (f1 − g1 )dτ + C1 ,
0
truncated after m = r , t
w2 (y, t) = (f2 + [w1 , f1 ]

r
m 0
ẏ = gm (y), (22) + [w1 , g1 ] − g2 )dτ + C2 ,
m! (28)
m=1 t
w3 (y, t) = f3 + [w1 , f2 ] + 2[w2 , f1 ]
the following theorem relates the solution of the r th 0
averaged system (22) and the complete averaged sys- + 2[w1 , g2 ]

tem (21). − [w1 , [w1 , g1 ]] + [w2 , g1 ] − g3 dτ
+ C3 ,
Theorem 2 (Theorem 6.5.2 in [18]) Let y(t) and y(t) and f1 = f(y, t, 0), f2 = 2 ∂f/∂|=0 , f3 =
be the solutions of (21) and (22), respectively. There 3 ∂ 2 f/∂ 2 |=0 , and f4 = 4 ∂ 3 f/∂ 3 |=0 .
exist positive constants K and 0 such that, for all  ≤
0
2.2 Higher-order averaging via chronological calculus
|y(t) − y(t)| = O( r ),
In this subsection, we present the higher-order averag-
for all 0 ≤ t ≤ K /. ing approach based on chronological calculus, which
was developed in the late 1970s by the Russian mathe-
maticians Agrachev and Gamkrelidze [28,39]. Chrono-
Therefore, we can express the solution of the original logical calculus is concerned with time-varying vec-
system (1) as tor fields and the asymptotic expansion generated by
their flows. Sarychev and Vela [27,31,34] utilized the
r −1 m
  chronological calculus tools to develop a generalized
x(t) = y(t) + wm (y(t), t) + O( r ). (23)
m! averaging theory (GAT), which gives arbitrarily high-
m=1
order approximation of the flow along a time-periodic
For the convenience of the reader, we report the first vector field. We provide below a brief description of
four terms of (22). the technique.
Consider the NLTP dynamics
1 T ẋ = f(x, t). (29)
g1 (y) = f1 dt, (24)
T 0
Equation (29) is similar to (1) with the only difference
1 T 
g2 (y) = f2 + [w1 , f1 ] + w1 , g1 dt, (25) that f does not depend explicitly on .2 The vector field
T 0 f is smooth in x and T -periodic in t. The solution of
1 T the differential Eq. (29) for the initial condition x(0) =
g3 (y) = f3 + [w1 , f2 ] + 2[w2 , f1 ] + 2[w1 , g2 ]
T 0 x0 is denoted by x(t) = φ(x0 , t) = φ f t (x0 ), where

+ [w1 , [g1 , w1 ]] + [w2 , g1 ] dt, (26) φ f
t : R → R is the t-dependent flow associated with
n n

the vector field f(x, t) and φ f

0 (x0 ) = φ(x0 , 0) = x0 .
1 T
g4 (y) = f4 + [w1 , f3 ] + 3[w2 , f2 ] Also, φ f
t is a diffeomorphism satisfying the differential
T 0 f
+ [w1 , [w1 , [w1 , g1 ]]] + equation φ̇ t = φ f
t ◦ f(x, t), where “◦” indicates the
composition of maps (functions).
− [w1 , [w2 , g1 ]] − 2[w2 , [w1 , g2 ]]
− 3[w1 , [w1 , g2 ]] 2 Although Sarychev [34] acknowledges the case of explicit
dependence of f on , it is the authors’ opinion that the anal-
+ [w3 , g1 ] + 3[w2 , g2 ] + 3[w1 , g3 ]
 ysis of such a case needs further investigation. Hence, we restrict
+ 3[w3 , f1 ]] dt, (27) our demonstration to systems in the form (29).

123
 M. Maggia et al.

We can write the solution of Eq. (29) using the To better understand the concept of GAT, we now
Volterra series expansion as recall some key results from the linear Floquet theory.
Consider the homogeneous LTP system,
x(t) = x0
∞
 t τ1 τm−1 ẋ = A(t)x, (32)
+ ··· Y(x0 ) dτm . . . dτ1
with A ∈ Rn×n and T -periodic in t. The flow asso-
m=1 0 0 0
(30) ciated with (32) can be represented by the so-called
 fundamental solution matrix, written as
where Y = Yτ1 ◦ Yτ2 ◦ · · · ◦ Yτm . Here, Yτ =  t 
∂Y →
f(x, τ ) and Y τ1 ◦ Y τ2 = ∂ xτ2 Y τ1 . The condi- Φ t = exp
A
A(τ )dτ ∈ Rn×n .
0
tions of convergence of series (30) are provided by
Agrachev and Gamkrelidze [29]. In a similar fash- The fundamental solution matrix of system (32), com-
ion to the exponential representation of solutions of puted after one period T , is known as the monodromy
time-invariant systems [30,40], Agrachev and Gamkre- matrix M = Φ A T.
lidze [29] expressed the flow of the time-varying vector Theorem 3 Linear Floquet Theorem [41, p. 11] Every
field f(x, t) as fundamental matrix solution Φ A
 t  t of the LTP system (32)
f → can be represented as the product Φ A t = P(t)e
Rt of a
φ t = exp f(x, τ ) dτ . (31) T -periodic matrix (i.e., P(t) = P(t + T )) and R is a
0
constant matrix given by R = T1 ln M.
In the chronological calculus, this flow is also known
as the right-chronological exponential. Theorem 4 System (32) is uniformly stable if and only
Then, the logarithm of (31) is defined as if all Floquet multipliers (eigenvalues of M) have mod-
uli less than 1.
Vt = ln(φ f
t ).

In other words, the flow along the time-invariant vector Based on Theorem 3, one can use the transformation
field Vt for a unit time is equivalent to the flow along x = Py to transform the LTP system (32) to the LTI
the time-varying vector field f(x, t) for a time t. system
It is noteworthy to mention that although Vt is an ẏ = Ry (33)
autonomous vector field, it is parametrized by time;
if the final time t is changed, the vector field will whose flow after one period is exactly equivalent to the
change. Furthermore, in Ref. [29] it was shown that flow of the original LTP system (32) after one period:
 (m) eRT = Φ A
Vt = ∞ m=1 V t , with T . That is, based on our discussion in the pre-
vious section on classical averaging, system (33) rep-
t τ1
Vt(m) = resents a stroboscopic averaged dynamics of the LTP
0 0 system (32). Moreover, one can determine the stabil-
τm−1
ity of the LTP system (32) from the Floquet multipliers
··· Gm (Yτ1 , . . . , Yτm )dτm . . . dτ1 ,
0 according to Theorem 4, or from the averaged dynamics
(33): The solution x(t) of (32) converges to the origin if
and Gm are the commutators polynomials. For instance,
and only if the solution y(t) of (33) converges to zero.
for m = 1, 2, 3, we obtain
Therefore, the LTP system (32) is stable if and only if
1 the averaged dynamics (33) is stable (i.e., the eigen-
G1 (ζ1 ) = ζ1 , G2 (ζ1 , ζ2 ) = [ζ2 , ζ1 ]
2 values of the R lie in the open left half of the complex
1 plane); there is no approximation here. It is complete
G3 (ζ1 , ζ2 , ζ3 ) =([ζ3 , [ζ2 , ζ1 ]] + [[ζ3 , ζ2 ], ζ1 ]) ,
6 averaging. The eigenvalues of R are called the Floquet
where [·, ·] denotes the commutator, which, for vec- exponents, or the characteristic exponents.
tor fields, coincides with the Lie bracket. There- Sarychev [38] utilized the Floquet concept of aver-
fore, chronological calculus provides an algorithmic aging without approximation to extend the Floquet
approach to determine the logarithm of time-periodic theory to the nonlinear case, equivalently extending
vector fields analytically. the averaging theorem to higher order. Similar to the

123
Reconciliation of two higher-order averaging techniques

Table 1 Linear versus nonlinear Floquet theory The vector field Λ in Eq. (34) can be written as a
Linear Nonlinear power series expansion in  such that
∞

ẋ = A(t)x ẋ = f(x, t) m
    ẋ = Λ(x) = Λm (x), (35)
→ → t m!
φ f
t
t = exp
ΦA 0 A(τ )dτ t = exp 0 f(x, τ ) dτ m=1
M = ΦA
T M = φTf where x is the average of x. According to Theorem 5,
R= 1
T ln M VT = 1
T ln M the solution of (29) oscillates around the solution of

ẏ = Ry ẋ = (x) = ∞ m=1
m
m! m (x)
the averaged NLTI system (35). Finally, if the series in
Eq. (35) is truncated after r terms, one arrives at the
r th-order averaged version

definition of the averaged dynamics R = T1 ln M = 

r
m
→ T ẏ = Λm (y). (36)
T ln Φ T = T ln exp 0 A(t)dt, Sarychev [38] intro-
1 A 1 m!
m=1
duced the notion of complete averaging of the NLTP
The power of this generalization is that the Λm
system (29) to indicate the autonomous vector field
terms can be computed analytically using Lie brack-
1 1 1 ets between the vector fields that describe the time-
Λ= VT = ln φ f
T = ln M, (34)
T T T periodic dynamics.
For the convenience of the reader, we report the first
where M = φ f T is the flow of f(x, t) computed after
four truncations as seen in Refs. [8,27,34].
one period T and is called the monodromy map, which
maps x0 to x(T ) (i.e., the solution after a period T ). The 1 T
Λ1 (y) = f dt, (37)
monodromy map can be seen as the nonlinear analo- T 0
 
gous to the monodromy matrix M in the linear Floquet 1 T t
Λ2 (y) = fτ dτ , f dt, (38)
decomposition. The vector field Λ is then equivalent to T 0 0
the matrix R. Furthermore, the linear Floquet theorem 3
Λ3 (y) = − T [Λ1 , Λ2 ]
can also be extended to the nonlinear case as follows 2
 t  t 
2 T
+ fτ dτ , fτ dτ, f dt
Theorem 5 Nonlinear Floquet Theorem (Theorem 3.2 T 0 0 0
  t  
in [34] and Theorem 8 in [27]) Assume that φ f t is 3 1 T 1 T
=− T f dt, fτ dτ , f dt
the flow generated by the NLTP system in (29). If the 2 T 0 T 0 0
 t  t 
diffeomorphism M = φ f T admits a logarithm VT , 2 T
+ fτ dτ , fτ dτ, f dt. (39)
then the flow φ f
t can be represented as a composition T 0 0 0
    
φ f
t = e Λt ◦ P of the flow eΛt of the autonomous vec-
t 2 T t τ τ1
Λ4 (y) = fτ2 dτ2 , fτ1 dτ1 , [fτ , f] dτ
tor field Λ = T1 VT and a T -periodic map P t (i.e., T 0 0 0 0
 t  τ  τ1   
P t = P t+T ).
+ fτ2 dτ2 , fτ1 dτ1 , fτ dτ, f
0 0 0
Moreover, the stability characteristics of this logarithm t τ  τ   
+ fτ1 dτ1 , fτ1 dτ1 , fτ , f dτ dt,
Λ are exactly related to the stability characteristics of 0 0 0
the original NLTP system (29) in a Floquet theorem (40)
fashion.
where f = f(y, t), fτ = f(y, τ ), fτ1 = f(y, τ1 ), and
Theorem 6 (Theorem 9 in [27]). If the monodromy fτ2 = f(y, τ2 ). It is worth mentioning that the expres-
map, M, of system (29) has a fixed point, then the flow, sions for Λ3 and Λ4 are not agreed upon between all
φ f
t , has a periodic orbit whose stability is determined cited references. In particular, there is a mismatch in
by the stability of the monodromy map. the sign of the term [Λ1 , Λ2 ] in Eq. (39) and in the sign
of the entire expression of Λ4 in Eq. (40). Based on
Theorem 6 clearly represents the nonlinear extension our analysis, we believe that the expressions for Λ3 in
of Theorem 4. Table 1 summarizes the parallel between Eq. (39) and for Λ4 in Eq. (40) are the correct ones. We
the linear and the nonlinear Floquet theorems. will provide further evidence in the next section.

123
 M. Maggia et al.

3 Reconciliation of the two higher-order averaging Lemma 1 Given two smooth vector fields X(x, t), Y
approaches (x, t) ∈ Rn , then
   
∂ ∂Y ∂X
In this section we will show that, under some key [X, Y] = X, + ,Y , (41)
∂t ∂t ∂t
assumptions, the two higher-order methodologies can
be proven equivalent despite the fact that they were and if X is time-invariant, the following are true
developed following very different and independent  
∂ ∂Y
derivations. We prove equivalence mathematically up [X, Y] = X, , (42)
∂t ∂t
to the fourth-order averaging (r = 4) under the follow-
ing assumptions: and
 
Assumption 1 Assume that the NLTP system is writ- [X, Y]dt = X, Ydt .
ten in the form (29)
ẋ = f(x, t), Proof Using Eq. (18), it follows that
with f smooth in x and T -periodic in t. ∂
[X, Y]
∂t
Assumption 2 (Stroboscopic averaging) Assume that  
∂ ∂Y ∂X
the arbitrary constants of integration in Eq. (17) and = X− Y ,
∂t ∂x ∂x
Algorithm 1 are taken to be null    
∂ ∂Y ∂Y ∂X ∂ ∂X ∂X ∂Y
= X+ − Y− ,
Cm = 0, for m ≥ 1. ∂t ∂x ∂x ∂t ∂t ∂x ∂x ∂t
   
∂ ∂Y ∂X ∂Y ∂Y ∂X ∂ ∂X
With Assumption 1, we reduce the analysis to the = X− + − Y,
∂x ∂t ∂x ∂t ∂x ∂t ∂x ∂t
class of systems in which the parameter  does not    
appear explicitly in f. As a result, Eq. (3) becomes ∂Y ∂X
= X, + ,Y .
f(x, t, ) = f(x, t) = f1 (x, t) (i.e., fm (x, t) = 0 for ∂t ∂t
m > 1) and system (1) is in the form (29). Other-
If X is time-invariant, then ∂X
∂t = 0 and Eq. (41) reduces
wise, the chronological calculus technique would not
to (42). Moreover,
be even applicable. The second assumption renders the  
first approach a fundamental concept similar to the ∂Y ∂X
[X, Y]dt = X− Y dt,
chronological calculus approach: The solution of the ∂x ∂x
averaged system coincides with that of the original sys- ∂Y ∂X
= Xdt − Ydt,
tem stroboscopically (after each cycle). This concept ∂x ∂x
was the basis of the chronological calculus averaging:  
∂ Ydt ∂X
The goal was to find an autonomous vector field whose = X− Ydt,
∂x ∂x
flow after one period matches that of the time-periodic  
vector field after one period. Therefore, while we pro- = X, Ydt .
vide below a mathematical proof for the equivalence
between the two approaches up to fourth order, this 
conceptual matching suggests that the equivalence can
be extended to all higher orders. Theorem 7 Given Assumptions 1 and 2, the following
First, we state the following properties of Lie are true: g1 = Λ1 , g2 = Λ2 , g3 = Λ3 , and g4 = Λ4 .
brackets which will be useful in reconciling the
Proof Given in the following Subsect. 3.1–3.4. 
two higher-order averaging techniques. Let X(x, t),
Y(x, t), Z(x, t) ∈ Rn and α ∈ R. The Lie brackets,
as defined in Eq. (18), satisfy the following properties: 3.1 First order
1. [X, Y] = −[Y, X] (skew-symmetry),
2. [X + Y, αZ] = α[X, Z] + α[Y, Z] (linearity), By comparing Eqs. (24) and (37), it trivially follows
3. [[X, Y], Z]+[[Y, Z], X]+[[Z, X], Y] = 0 (Jacobi that g1 = Λ1 . This equivalence holds also without
identity). Assumption 2.

123
Reconciliation of two higher-order averaging techniques

3.2 Second order 3.3 Third order

For the second-order terms, we start by expanding g2 In a similar fashion to the second-order averaging, let
given in Eq. (25) as us expand g3 by substituting the values of w1 and w2
  t   t 
1 T found in Eq. (28), into Eq. (26).
g2 = f2 + f1 dτ, f1 − g1 dτ , f1  
  t 
T 0 0 0 1 T t
 t   t  g3 = f3 +
f1 dτ, f2 − g1 dτ, f2
T 0 0 0
+ f1 dτ , g1 − g1 dτ , g1 (43)  t   t 
0
 0 +2 f2 dτ, f1 + f2 dτ, g1
0 0
+ [C1 , f1 ] + [C1 , g1 ] dt.  t τ  
+2 f1 dτ1 , f1 dτ, f1
From Assumption 1, we have f2 = 0 and f1 = f, 0
 t τ
0
 
and from Assumption 2, C1 = 0. Equation (43) then + f1 dτ1 , f1 dτ, g1 +
becomes
 t   t  0
 t τ
0
 
1 T
g2 = f dτ, f − g1 dτ , f −2 g1 dτ1 , f1 dτ, f1
T 0 0 0 0 0
 t   t   t τ  
+ f dτ , g1 − g1 dτ , g1 dt. − g1 dτ1 , f1 dτ, g1
0 0
0 0  t τ  
Since g is time-invariant, Lemma 1 implies +2 f1 dτ1 , g1 dτ, f1
 t 1 
  0 0
g1 dτ , g1 = t g1 , g1 = t g1 , g1 = 0.  t τ  
0 + f1 dτ1 , g1 dτ, g1
To recover the second-order averaging using chrono- 0
 
0

t T
logical calculus, Λ2 in Eq. (38), it remains to show that −2
1
 t   t  T
f2
1 T 
0 0
  
f dτ , g1 − g1 dτ , f dt = 0. t
T 0 0 0 + f1 dτ, f1 dt dτ, f1
(44) 
0
t 1 T 
Using the skew-symmetry property of the Lie bracket − f2
(18), the left-hand side of Eq. (44) becomes 0 T 0
 t     t   
1 T t
+ f1 dτ, f1 dt dτ, g1
f dτ, g1 + f, g1 dτ, dt. 0t  t 
T 0 0 0
Applying property (41), we obtain − f1 dτ, f1 dτ, g1
 t  0
 t
0
 t 
1 T ∂ t
f dτ, g1 dτ, dt. + g1 dτ, f1 dτ, g1
T 0 ∂t 0 0 0 0
Using the fundamental theorem of calculus, the expres-  t T   t  
1
sion above simplifies to +2 f1 dτ, f2 + f1 dτ, f1 dt
T
 T T  
0
 
0
 
0
1 t 1 t T
f dt, g1 dt, . −2 f2 +
g1 dτ, f1 dτ, f1 dt
T 0 0 0 0 T 0
T T  t   t 
Finally, noting that 0 f dt = T g1 , and 0 g1 dt =
+2 [C1 , f1 ] dτ, f1 − 2 [C1 , g1 ] dτ, f1
T g1 , 0 0
 t 
1 
T g1 , T g1 = 0. + [C1 , f2 ] + (f1 − g1 ) dτ + 3[C1 , g1 ]
T 0
Hence, the second-order averaging via Lie transform  t 
given in Eq. (43) reduces to + [C1 , g1 ] dτ , g1
 t  
0
 t 
1 T
g2 = f dτ , f dt. (45) + C1 , − (f1 − g1 ) dτ, g1
T 0 0 0
Comparing Eqs. (45) and (38), it is easy to see that  T   t  
1
g2 = Λ2 . + 2 C1 , f2 + f1 dτ , f1 dt
T 0 0

123
 M. Maggia et al.

  t  t 
+ [C1 , 3[C1 , g1 ]] dt + 3[C2 , g1 ]. (46) 2 T
+ f dτ , f dτ, f dt.
T 0 0 0
According to Assumption 1, f = f1 and f2 = f3 = In order to prove that g3 = Λ3 , we now need to show
0. Assumption 2 leads to C1 = C2 = 0. Therefore, that
Eq. (46) simplifies to  t τ  
T  2 T
 t τ   g1 , g2 − f, f dτ1 , f dτ dt =
2 T 2 T 0
g3 = f dτ1 , f dτ, f dt 0
 t
0
 t 
T 0 0 0 3T  2 T
T  t τ   − g1 , g2 + f dτ, f dτ, f dt,
2 2 T 0
− g1 dτ1 , f dτ, f dt 0 0
T 0 0 0 or equivalently, that
 t τ   
1 T  2 T t τ  
+ f dτ1 , f dτ, g1 dt 2T g1 , g2 = f, f dτ1 , f dτ dt
T 0 0 0 T 0
 t τ   
0 0
 t  (49)
1 T 2 T t
− g1 dτ1 , f dτ, g1 dt + f dτ, f dτ, f dt.
T 0 0 0 T 0
T  t τ   0 0
2 Let us consider the right-hand side of Eq. (49)
+ f dτ1 , g1 dτ, f dt 
T 0 0 0 2 T t τ  
 t τ   f, f dτ1 , f dτ dt
1 T T 0
+ f dτ1 , g1 dτ, g1 dt 0
 t
0
 t 
T 0 0 0 2 T
  (47) + f dτ, f dτ, f dt.
2 T t
T 0
− g2 dτ, f dt 0 0
T 0 0 Using the linearity of integration and the Lie bracket
 t 
1 T property (41), we obtain
−  t t τ  
g2 dτ, g1 dt
T 0 0 2 T ∂
T  t  t  f dτ, f dτ1 , f dτ dt,
1 T 0 ∂t 0
− f dτ, f dτ, g1 dt 0 0
T 0 0 0 which is further simplified, by applying the fundamen-
 t  t 
1 T tal theorem of calculus, to
+ g1 dτ, f dτ, g1 dt  T T  t  
T 0 0 0 2
 t  f dτ, f dτ, f dt .
2 T T 0 0 0
+ f dτ, g2 dt
T 0 0 Finally, recalling the expressions for g1 and g2 in
T  t  Eqs. (24) and (45), respectively, we obtain
2 
− g1 dτ, g2 dt,
T 0 0 2T g1 , g2 .
where g2 is given by Eq. (45). Using the Lie brackets Thus, g3 = Λ3 and the equivalence of the third-order
properties and some other algebraic manipulations,3 terms of the two higher-order averaging techniques is
Eq. (47) can be further simplified to proved.
T 
g3 = g1 , g2
2
 t τ   3.4 Fourth order
2 T
− f, f dτ1 , f dτ dt. (48)
T 0 0 0 According to Assumption 1, f = f1 and f2 = f3 = f4 =
Using the fact that g1 = Λ1 and g2 = Λ2 , we rewrite 0. Assumption 2 leads to C1 = C2 = C3 = 0. Sub-
Eq. (39) as, stituting the expressions for w1 , w2 and w3 in Eq. (28)
3T  into Eq. (27) and after some algebraic manipulations,
Λ3 = − g1 , g2 the expression for g4 in Eq. (27) reduces to
2
T  t  τ  τ2   
3 6
The process of reducing Eqs. (47)–(48) is rather lengthy and g4 = f dτ2 , f dτ1 , f dτ, f dt
would take several pages. Thus, the authors decided not to include T 0 0 0 0
(50)
this part in this paper and refer the reader to our Mathematica® T2  
file [42] for the complete proof. − g1 , g1 , g2 + T g1 , g3 ,
2

123
Reconciliation of two higher-order averaging techniques

where g1 , g2 and g3 are given by Eqs. (24), (45) 4 Applications

and (48), respectively. 
It can be shown that g4 = Λ4 , or equivalently, In this section, we present two applications of higher-
order averaging to mechanical systems. First, we apply
    higher-order averaging to a classical problem, namely
6 T t τ
the Kapitza inverted pendulum problem (see, for exam-
hτ1 dτ1 , f dτ, f dt
T 0 0 0 ple, [43]). Then we consider a modern application
T 2   of flapping flight of biological flyers and their man-
− g1 , g1 , g2 + T g1 , g3 made counterparts: flapping-wing micro-air-vehicles
2
 t τ  (see [22,44]).
2 T
= hτ1 dτ1 , [fτ , f] dτ (51)
T 0 0 0
 t τ  
4.1 Kapitza’s inverted pendulum
+ hτ1 dτ1 , fτ dτ, f
0 0
t τ   The Kapitza inverted pendulum problem has been stud-
+ fτ1 dτ1 , [hτ , f] dτ dt, ied for decades because it gives interesting insights into
0 0
the concept of vibrational stabilization. In fact, it has
been shown that high-frequency, high-amplitude, peri-
τ  τ  odic forcing, applied in the vertical direction to the pen-
where hτ1 = 0 1 fτ2 dτ1 , fτ1 , and hτ = 0 fτ1 dτ, fτ .
dulum pivot, stabilizes the pendulum about its inverted
For the details of the proof of the equivalence in Eq. (51)
position [10,20]. This problem was first analyzed by
we refer once again to our Mathematica® file [42].
Stephenson [45] in 1908 using a periodic series solution
With this reconciliation, each approach can bene-
to the linearized system (Hill’s equation) and then by
fit from the tools and results developed in the other
Kapitza [43,46] considering the full nonlinear dynam-
approach. For example, there are no recursive formulas
ics using an averaged potential approach. The reader
for higher-order averaging using chronological calcu-
may like to watch the interesting video in Ref. [47]
lus: Neither Vela [27] nor Sarychev [34] provided such
that clearly shows the phenomenon and its associated
recursive formulas. For a specific desired r th-order
stabilizing stiffness mechanism.
averaging, one has to derive it following the approaches
The pendulum dynamics, linearized about its unsta-
presented by Sarychev [34], Vela [27], and the one dis-
ble vertical equilibrium, can be written as
cussed above. Instead, Algorithm 1 provides a straight-
forward recursive formula for arbitrary higher-order θ̈ + (−Ω + B1 cos(ωt))θ = 0, (52)
averaging terms using Lie transform. As such, with the where θ and θ̇ represent the angular position and veloc-
proved equivalence, one can use Algorithm 1, setting ity, respectively, ω is the forcing frequency, B1 is the
f1 = f, fm = 0 for m > 1 and Cm = 0 for all m, forcing amplitude, and Ω = g/l (g is the gravitational
to determine a recursive formula for the higher-order acceleration and l the pendulum length). Figure 2 shows
averaging terms defined through chronological calcu- a schematic of the inverted pendulum oscillating ver-
lus. tically about its pivot. Setting x1 = θ and x2 = θ̇,
On the other hand, Sarychev [34] and Vela [27] Eq. (52) can be written in state space representation as
showed that system (35) represents complete averaging    
ẋ1 x2
of the NLTP system (29). In other words, the full series = .
ẋ2 (Ω − B1 cos(ωt))x1
in Eq. (35) represents the logarithm of the Monodromy
The problem of interest (vibrational control problem) is
map. Therefore, if this series converges to Λ, the sta-
for high amplitude and high frequency; thus, we intro-
bility characteristics of this logarithm Λ are exactly
duce the scaling B1 = B/ and ω = 1/, where  is a
related to the stability characteristics of the original
small parameter. This yields
NLTP system (29). Through the proved equivalence,    
one can apply Theorem 6 to the full series (21) of the ẋ1 x2
= . (53)
averaged dynamics using Lie transform. This result is ẋ2 (Ω − B cos( t ))x1
important because it has not been proved within the System (53) is not amenable to the averaging theorem: f
classical averaging formulation. is not analytic in epsilon because of the high-amplitude

123
 M. Maggia et al.

x(t)

g f +g g
l

θ F y(0) = z(t)

x(0) = z(0) = x0
B1 cos(ωt)
Fig. 3 The flow along f + g is equivalent to the flow along the
vector fields F and then g

and the first-order averaged dynamics of the entire sys-

Fig. 2 Kapitza’s inverted pendulum tem would be obtained by averaging F.
To apply the VOC formula to system (53) with
f(x, t) = (x2 , Ω x1 )T and g(x, t) = (0, −(B/)
forcing term proportional to 1/. This scenario is typ-
cos(t/)x1 )T , we first obtain the flow of g
ical with vibrational control systems (high-amplitude  
periodically forced systems). The standard remedy for g z1 
φ t (z) =
this issue is performed by applying the variation of con- z 2 − B sin t z 1
stants (VOC) formula [7].   
1 0 z1
=
−B sin( t ) 1 z2
4.1.1 The nonlinear variation of constants formula
then
 
The VOC formula is used to decouple the flow along g z 2 − B sin( t )z 1
f ◦ (φ t )(z) = ,
two vector fields. In particular, when applied to the Ωz 1
multi-scale system (53), it decouples it into two sys- and since
tems, each of which is amenable to the averaging theo-  g −1  
∂φ t 1 0
rem. The VOC is quite general; it is applicable to sys- = ,
∂z B sin( t ) 1
tems in the form
we determine the pullback vector field F as
ẋ = f(x, t) + g (x, t) . (54)   
1 0 z1
The VOC formula decouples system (54) into the two F(z, t) =
B sin( t ) 1 z 2 − B sin( t )z 1
systems   (56)
z 2 − B sin(

t
 )z 1
ż = F(z), z(0) = x0 , = .
(55) Ω − B 2 sin2 ( t ) z 1 + B sin( t )z 2
ẏ = g(y, t), y(0) = z(t),
Finally, we perform the scaling τ = ωt to system (56),
where F is the pullback of f along the flow of g, defined to render it in the averaging-canonical form (weakly
as forced system with f analytic in ).
 g −1 
g ∗  ∂φ t g z = F(z, τ )
F(z, t) = (φ t ) f (z) = ◦ f ◦ φ t (z),  
∂z z 2 − B sin(τ )z 1 (57)
=  ,
g Ω − B 2 sin2 (τ ) z 1 + B sin(τ )z 2
and φ t is the flow along the vector field g for time t.
Figure 3 shows the idea behind the VOC formula: The where z = dz/dτ .
flow along f +g corresponds to the flow along the vector We apply higher-order averaging to system (57). In
field F, and then g. The two systems in Eq. (55) should order to assess the accuracy of the higher-order averag-
be amenable to the averaging theorem, perhaps after ing approximations, we compare the obtained results
a simple transformation. Moreover, since g is a zero- with the Floquet stability map, obtained by numer-
mean, time-periodic vector field, its averaging vanishes ical application of the Floquet theorem to the LTP

123
Reconciliation of two higher-order averaging techniques

Fig. 4 Stability map for the 1

Kapitza pendulum. The 1st & 2nd
shaded region represents the 0.9 3rd
Floquet stable region. The 4th
colored curves represent the 0.8 5th
higher-order stroboscopic
averaging approximations 0.7
(with VOC formula) of the
lower stability boundary 0.6

0.5

0.4

0.3

0.2

0.1

0
-0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0

system (53). The approximations are computed up to as the frequency is increased, which matches the intu-
fifth order using the algorithm described in Sect. 2.1, ition behind averaging. These two conclusions indi-
and by choosing the stroboscopic averaging, namely cate that, for a required level of accuracy, the thresh-
in Algorithm 1, we set Ck = 0 for k = 1, 2, 3, 4. old frequency decreases with the averaging order. That
Hence, the two averaging techniques will lead to iden- is, higher-order averaging allows stabilization at fre-
tical results due to the proved equivalence in this case quencies lower than that demanded by direct averaging,
of ∂F/∂ = 0. Figure 4 illustrates the stability map hence making vibrational control more feasible.
of the Kapitza pendulum in a similar fashion to the Figure 6 shows the angular displacement and veloc-
work of Berg and Wickramasinghe [20]: The plotted ity (θ , θ̇) of the Kapitza pendulum with initial con-
quantities are defined as α = −Ω/ω2 and β = B/ω, ditions θ0 = 30◦ and θ̇0 = 0◦ /s and for Ω = 2,
and the region colored in yellow represents the Floquet B = 3, and ω = 20 rad/s. The three curves repre-
stable region. In addition, the figure shows the approx- sent the integration of the complete LTP dynamics (57),
imations of the lower stability boundary obtained via the first-order and third-order averaged dynamics. The
higher-order averaging. Although only the lower sta- solutions of the averaged dynamics have been trans-
bility boundary is captured, it is evident that the accu- formed according to Eq. (23). While both the averaged
racy of the approximations increases with the averag- solutions correctly capture the stability of the origi-
ing order. That is, the higher-order averaging stabil- nal system, the third-order approximation is practically
ity boundary converges asymptotically to the Floquet indistinguishable from the solution of the full LTP sys-
boundary as the averaging order increases. tem.
Figure 5 shows the percent error between the Floquet We also present in Figs. 7 and 8 the results for a non-
lower stability boundary and the higher-order averag- stroboscopic averaging case. In particular, the integra-
ing approximations as the forcing frequency ω is var- tion constants Ck are computed using Eq. (20). Similar
ied. Two important conclusions can be drawn from this conclusions to the stroboscopic case can be drawn.
figure. First, for a fixed ω, the error decreases as the
averaging order is increased, which matches the intu-
ition behind the asymptotically convergent series as 4.1.2 Higher-order averaging without VOC
more terms are taken [15]. Second, the error associated
with all the averaging order approximations decreases As shown above, vibrational control systems are typ-
ically forced by high-amplitude periodic inputs and

123
 M. Maggia et al.

Fig. 5 Percent error in the 0.8

lower stability boundary for 1st & 2nd
Ω = 2 for the Kapitza 3rd
pendulum using 0.7 4th
higher-order stroboscopic 5th
averaging with the VOC 0.6
formula

0.5

Error (%)
0.4

0.3

0.2

0.1

0
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

30
the time-periodic terms, which may not be analytically
20
tractable for many systems. Here we show that higher-
10 order averaging alone (without the VOC) may capture
0 some of the dynamical features that are lost by disre-
-10 garding the VOC. In order to demonstrate this point,
-20 consider again the system in Eq. (53) after the transfor-
mation τ = ωt = t/
-30
0 5 10 15 20 25 30    
x1 x2
= , (58)
x2 (Ω − ωB cos(τ ))x1
50
Figure 9 shows the percent error in the stability
boundaries obtained by applying higher-order averag-
ing directly (without VOC) to system (58). Similar con-
0 clusions to those of Fig. 5 (with the VOC formula) can
be drawn: The accuracy increases as ω increases (i.e.,
 decreases) and/or the averaging order is increased.
-50 However, in comparison with the higher-order averag-
0 5 10 15 20 25 30
ing with VOC, the accuracy of the approximations in
Fig. 9 is noticeably lower; a significantly higher-order
Fig. 6 Integration of the full LTP dynamics, first- and third-order averaging may be required to attain an acceptable level
averaged dynamics of the Kapitza pendulum, after applying the of accuracy. Finally, Fig. 10 shows the angular displace-
VOC formula. The simulations are obtained for Ω = 2, B = 3, ment and velocity of the Kapitza pendulum with initial
ω = 20 rad/s, and initial conditions θ0 = 30◦ and θ̇0 = 0◦ /s
conditions θ0 = 30◦ and θ̇0 = 0◦ /s and for Ω = 2,
B = 3, and ω = 20 rad/s. The four curves represent
the integration of the complete LTP dynamics (58), the
therefore are not in the averaging-canonical form; the first-order, third-order and fifth-order averaged dynam-
VOC formula is necessary to ensure a rigorous appli- ics. The solutions of the averaged dynamics have been
cation of the averaging theorem. However, as shown adjusted according to Eq. (23). In this case, the solution
above, the VOC requires computing the flow map along of the first-order averaged dynamics provides wrong

123
Reconciliation of two higher-order averaging techniques

Fig. 7 Stability map for the

Kapitza pendulum. The 1st & 2nd
shaded region represents the 0.9 3rd & 4th
5th & 6th
Floquet stable region. The
0.8 7th & 8th
colored curves represent the
higher-order
non-stroboscopic averaging 0.7
approximations (with VOC
formula) of the lower 0.6
stability boundary
0.5

0.4

0.3

0.2

0.1

-0.3 -0.25 -0.2 -0.15 -0.1 -0.05

Fig. 8 Percent error in the

lower stability boundary for 1st & 2nd
Ω = 2 for the Kapitza 3rd & 4th
5th & 6th
pendulum using 0.125
7th & 8th
higher-order
non-stroboscopic averaging
with the VOC formula
0.1
Error (%)

0.075

0.05

0.025

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

information about the stability of the original system, when the application of the VOC formula is not possi-
while the stability characteristics are correctly captured ble, such as in the modern example shown below.
by the higher-order averaged dynamics. As expected,
when the VOC formula is not applied, a higher aver-
4.2 Hovering flight of insects and flapping-wing
aging order (fifth order versus third order for the VOC
micro-air-vehicles (FWMAVs)
case) may be required to obtain a close approximation
of the solution of the original LTP system.
Flapping flight dynamics is a very rich dynamical sys-
In conclusion, higher-order averaging might be ben-
tem because it is typically represented by a multi-body,
eficial to the stability analysis of complex systems even
nonlinear, time-varying model. To simplify the analy-

123
 M. Maggia et al.

Fig. 9 Percent error in the 14

lower stability boundary for 3rd & 4th
Ω = 2 for the Kapitza 5th & 6th
pendulum using 12 7th & 8th
higher-order averaging 9th & 10th
without the VOC formula
10

Error (%)
6

0
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

150
flapping flight dynamics in the longitudinal plane (x −z
plane in Fig. 11) will be governed by the exact same
100
set of equations governing conventional airplanes [49]
⎛ ⎞
50
⎛ ⎞ ⎛ ⎞ 1
X (x, t)
u̇ −qw − g sin θ ⎜ m ⎟
0 ⎜ ẇ ⎟ ⎜ qu + g cos θ ⎟ ⎜ m1 Z (x, t) ⎟
⎜ ⎟=⎜ ⎟+⎜ ⎟
⎝ q̇ ⎠ ⎝ 0 ⎠ ⎜ 1 M(x, t) ⎟ , (59)
-50 ⎝ Iy ⎠
0 1 2 3 4 5 θ̇ q 0
200
where I y and m are the pitch inertia and the body
100 mass, respectively, and g is the gravitational accel-
eration. Figure 11 illustrates a flapping-wing micro-
0 air-vehicle (FWMAV) in the longitudinal plane (x − z
plane) where the state vector x includes the forward and
-100
normal velocity components u, w, the pitching angular
velocity q, and the pitching angle θ . In Eq. (59), X , Z
-200
0 1 2 3 4 5 are the forward and normal aerodynamic forces, respec-
tively, and M is the aerodynamic pitching moment.
Fig. 10 Integration of the full LTP dynamics, first-, third- and Details of the aerodynamic model (dependence of
fifth-order averaged dynamics of the Kapitza pendulum, without the aerodynamic forces X , Z and M on the state vari-
applying the VOC formula. The simulations are obtained for ables u, w and q) are taken from Ref. [50] and given in
Ω = 2, B = 3, ω = 20 rad/s, and initial conditions θ0 = 30◦
and θ̇0 = 0◦ /s
“Appendix A.”
The main distinction between conventional air-
planes and FWMAVs or insect flight is that the aero-
sis of such a complicated system, it may be reasonable dynamic forces (X , Z and M) are essentially time-
to neglect the wing flexibility and inertial effects, even varying in the latter case, which renders system (59)
though it is a controversial assumption among the flap- time-varying nature. That is, system (59) can be writ-
ping flight dynamics community [6,48]. Doing so, the ten as

123
Reconciliation of two higher-order averaging techniques

ping flight dynamics of insects and micro-air-vehicles

[52–61]. Using the crude averaging approach, it was
shown that the pitch dynamics of hovering insects and
FWMAVs are unstable due to lack of pitching stiffness
[54–59]. However, we recall the expressive quote from
Sanders and Verhulst book [17] presented in Introduc-
tion section stressing that averaging has to be rigorously
justified. In fact, the flapping flight dynamical system
as a high-amplitude periodically forced system is not
directly amenable to the averaging theorem. Moreover,
its dynamics are too complicated to allow analytical
Fig. 11 Longitudinal view of a hovering FWMAV application of the VOC. Therefore, there is a need for
higher-order averaging [44].
Prof. Nayfeh had a long career studying energy
ẋ = f(x) + ga (x, t), (60)
transfer between high-frequency modes and low-
where the autonomous vector field f represents the iner- frequency modes [62–67]. It is noteworthy to mention
tial and gravitational loads and the non-autonomous that Prof. Nayfeh opposed direct (crude) averaging of
vector field ga represents the aerodynamic loads. In flapping flight dynamics despite the seemingly large
conclusion, unlike conventional airplanes, flapping separation between the system’s two timescales that
flight dynamics (i.e., (59) and (60)) is represented by intuitively invokes averaging. Hence, this whole line
a nonlinear, time-periodic (NLTP) system whose sta- of research of higher-order averaging of insect flight
bility can be studied using one of the two approaches dynamics was initiated. In an earlier work [68], the
shown in Fig. 1. third author and Nayfeh utilized the multiple scales
Hovering is mostly associated with relatively higher method [15,32] to analyze the hovering flight dynamics
flapping frequencies in comparison with forward flight. of insects and FWMAVs up to second-order accuracy;
The flapping frequencies of hovering insects are usually they indeed highlighted the shortcomings of crude aver-
in the 20–1000 Hz range [51]. Therefore, the dynam- aging and showed that the parametric excitation due to
ics of hovering insects exhibit two timescales. The fast the oscillatory aerodynamic loads may naturally (with-
timescale is associated with the flapping motion and out feedback) stabilize the insect flight dynamics.
the aerodynamic loads, while the slow one is associ- In this subsection, we present how high-order aver-
ated with the aggregate motion of the body. The contrast aging captures this important dynamical feature in the
between the fast and slow timescales can be clearly seen hovering flight dynamics of insects and FWMAVs [44].
if one tries to observe a flying insect: A bare human eye In this demonstration, we adopt the flight dynamic
can track the trajectory of an insect body in the space, model (i.e., model for the aerodynamic vector field ga
but can hardly track the motion of its wings, which in terms of the state variables) developed earlier in [50].
move much faster. Therefore, we definitely have two Strictly speaking, system (60) (equivalently (59)) is
timescales. A natural question is to ask about the ratio not in the averaging-canonical form because the time-
between these two timescales. One of the slowest flap- periodic aerodynamic vector field ga (x, t) includes
ping insects, the hawkmoth, has its flapping frequency high-amplitude terms. Therefore, a rigorous applica-
(i.e., forcing frequency) about 30 times its natural fre- tion of the averaging theorem necessitates application
quency of flight dynamics (specifically the short-period of the VOC formula before averaging. However, it is
flight dynamics) [50,51]. This large separation between not possible to analytically compute the flow along
the two timescales invokes averaging; the body of an the NLTP vector field ga . This precludes application
insect hovering over a flower seems to oscillate in all of direct averaging, which is not clear until a care-
directions. Yet, on average, the insect is hovering over ful scaling is assigned to each term in ga . Yet, direct
the flower. Therefore, averaging seems a very intuitive (crude) averaging has been ubiquitously applied to ana-
tool to analyze flapping flight dynamics. lyze flapping flight dynamics over the past two decades,
Over the last two decades, direct (crude) averaging relying on the seemingly large separation between the
has been the standard analysis technique of the flap- two timescales [52–61].

123
 M. Maggia et al.

Table 2 Ratios of flapping frequency to natural frequency ω/ωn and eigenvalues assessing stability (using first- and second-order
averaging) for five insects
ω
Insect ωn λ1st λ2nd

Hawkmoth 28.78 [−11.89, − 3.30, 0.19 ± 5.74i] [−10.40, − 3.09, − 0.66 ± 3.72i]
Crane fly 50.62 [−47.71, − 17.31, − 1.13 ± 5.53i] [−45.76, − 16.26, − 13.16, 7.90]
Bumblebee 144.46 [−11.63, − 4.39, 1.58 ± 6.55i] [−11.26, − 4.37, 1.38 ± 6.17i]
Dragonfly 145.50 [−13.11, − 7.03, 1.34 ± 6.65i] [−12.56, − 6.98, 1.04 ± 5.99i]
Hoverfly 113.98 [−14.01, − 7.27, 2.13 ± 8.56i] [−13.37, − 7.24, 1.79 ± 7.92i]

As explained by Khan and Agrawal [60], the time systems are qualitatively different. Moreover, we can
variable in Eq. (60) is scaled as τ = ωωn t, where ωn is the even identify the vibrational stabilization mechanism
natural frequency of the body flight dynamics (short- comparing the matrices of the linearized direct (first-
period mode) and ω is the flapping (forcing) frequency. order) and second-order averaged systems for one of
When ωωn is large enough and assuming that all the the insects (e.g., hawkmoth)
terms in f and ga are equally scaled (which is not true in
general), the dynamics expressed in terms of τ are in the ⎡ ⎤
−3.59 0 0 −9.81
averaging-canonical form with  = ωωn . It is important ⎢ 0 −3.30 0 0 ⎥
D(Λ1 )(0) = ⎢
⎣ 39.95
⎥,
to note that, for hovering insects with relatively low 0 −7.92 0 ⎦
flapping frequencies (e.g., for the hawkmoth, ωωn ≈ 30) 0 0 1 0
⎡ ⎤
the ratio ωωn has usually been considered high enough to −3.58 0 0 −9.81
⎢ 0 −3.09 0 0 ⎥
justify direct averaging (see, for example, [50,52,55]). D(Λ1 +  2 Λ2 )(0) = ⎢
⎣ 29.98
⎥.
0 −8.13 −28.45 ⎦
The averaged dynamics of Eq. (60) are written (see
−2.90 0 0.96 0
[50]) as
ẋ = f(x) + ga (x) (61)
As can be seen, the element (3, 4), which corresponds
where x is the averaged state vector and ga (x) repre- to pitching acceleration (or pitching moment) due to
sents the ga (x, t) averaged aerodynamic loads: ga (x) = a pitching angle, is changed from 0, in the crude-

1 T
T 0 ga (x, τ ) dτ . averaged system, to −28.45 in the second-order aver-
After applying second-order averaging to system aged system. That is, a significant pitch stiffness (whose
(60), we obtain a significant change in the stability lack was the main reason behind the claimed instabil-
characteristics in comparison with crude averaging. ity) is created due to high-frequency, high-amplitude,
Table 2 shows the ratios ωωn for five different insects and zero-mean terms in ga which is a typical vibrational
the eigenvalues λ1st and λ2nd , computed after apply- stabilization behavior. Figure 12 shows a simulation
ing first- and second-order averaging, respectively. The of the hawkmoth dynamics for an initial θ disturbance
eigenvalues are used to assess stability. of 30◦ . The curves represent the simulation of the full
Some important observations can be made about NLTP system, the first-order averaged dynamics and
the results in Table 2. First, for very-high-frequency the second-order averaged dynamics. The solution of
cases, first-order (direct averaging) and second-order the second-order averaged dynamics greatly improves
averaging give the same information about the aver- the approximation obtained by applying direct averag-
aged system and its stability such as in the cases of the ing only.
bumblebee, the dragonfly and the hoverfly. In fact, for This example shows the power of higher-order aver-
these insects the ratio ωωn is very large (larger than 100). aging in capturing hidden stabilization mechanisms in
Also, this result is consistent with that obtained by the dynamically rich systems in an analytic manner that
method of multiple scales [68]. On the other hand, when allows deep study of the underpinning physics. This
the ratio ωωn is relatively lower, such as in the case of scrutiny of the physics behind vibrational stabilization
the crane fly and of the hawkmoth, the stability char- in flapping flight and its experimental validation using
acteristics of the first-order and second-order averaged motion data from a real hawkmoth is presented in [69].

123
Reconciliation of two higher-order averaging techniques

1 0.4
to a natural stabilizing action without feedback. In addi-
0.5 0.2 tion to these applications which demonstrate the power

w (m/sec)
u (m/sec)

0 0 of higher-order averaging, we show the equivalence

-0.5 -0.2
between two higher-order averaging approaches that
were developed independently in the literature within
-1 -0.4
0 50 Full NLTP
100 Dynamics0 50 100 two different communities and using different analysis
Second-Order Averaged Dynamics
First-Order Averaged Dynamics tools: the classical averaging approach using pertur-
4 40
bation theory and a more modern averaging approach
2 20
using chronological calculus. We show how the use of
q (rad/sec)

0 0 chronological calculus and Lie algebraic tools leads to

-2 -20 extension of the Floquet theorem to nonlinear systems,
-4 -40
equivalently extending the averaging theorem to arbi-
0 50 100 0 50 100 trarily high orders. We prove equivalence between the
t/T t/T
two averaging approaches mathematically up to fourth
Fig. 12 Integration of the full NLTP dynamics, first- and second- order and also show that the two approaches are con-
order averaged dynamics for the hawkmoth. The simulations are ceptually equivalent for all orders: Both rely on the
obtained for an initial disturbance of θ0 = 30◦ fact that the flows of the averaged system and the time-
periodic system match after a complete cycle. Based
on this equivalence, we conclude that the full series
5 Concluding remarks solution of the classical averaging approach represents
complete averaging in the sense that its exponential sta-
In this paper, we review the two main techniques to bility implies exponential stability for the original time-
analyze nonlinear time-periodic systems: the numer- periodic system for all values of the “small” parame-
ical approach using Floquet theorem and the analyt- ter (i.e., no large separation between the timescales is
ical approach using the averaging theorem with par- needed). Such a result was not developed in the classi-
ticular focus on the second approach. We discuss the cal averaging literature. On the other hand, the recursive
technical issues with the classical averaging theorem: formula developed in the classical averaging literature
the requirement of a weakly forced system and/or a can be directly used to determine higher-order aver-
large separation between the system’s two timescales. aging terms in the chronological calculus approach,
These requirements confine the applicability of the neat which was not presented in its respective literature. This
averaging approach to a small class of systems that article is considered as a comprehensive introductory
excludes vibrational control systems, which are typi- reference on the analysis of time-periodic systems that
cally high-amplitude periodically forced systems. We reviews and unifies the classical averaging theorem, the
present higher-order averaging as a remedy for these Floquet theorem, the higher-order averaging technique
issues. It is found that higher-order averaging captures using perturbation theory, and that uses chronological
more dynamical features than direct (first-order) aver- calculus and Lie algebraic tools. We show connections,
aging. In particular, we find that the error in the stabil- equivalence, similarity and differences between these
ity boundary decreases asymptotically as the averaging approaches.
order increases. We demonstrate this behavior on the
classical example of the Kapitza pendulum: inverted Acknowledgements This work is dedicated to the memory
of Prof. Ali Hassan Nayfeh whose contributions in the field
pendulum whose pivot is subject to vertical oscillation. of nonlinear dynamics and averaging of time-periodic systems
We also present a modern example with richer dynam- are immense. His passion and engineering intuition were always
ics: flapping flight dynamics of insects and micro-air- remarkable and inspiring. To honor his memory, we dedicate this
vehicles. Over the past two decades, direct averaging effort, which he enlightened its path before leaving. The authors
are thankful to the fruitful discussions with Dr. Muhammad Hajj,
has led to the conclusion that insects are unstable at Dr. Craig Woolsey, Dr. Jordan Berg, Dr. Patricio Vela, Dr. Sevak
hover. In contrast, we show that higher-order averaging Tahmasian, Dr. Ahmed Hassan and the support of the National
indicates vibrational stabilization in the form of a pitch
stiffness. That is, the high-frequency, zero-mean forces,
which are typically neglected by direct averaging, lead

123
 M. Maggia et al.

Science Foundation Grant CMMI-1709746 and its continuation

CMMI-1846308. 2Δx
Mq = − |ϕ̇| cos ϕ cos η (K 12 x h + K 22 sin ϕ)
Iy
Compliance with ethical standards 1 
+ ϕ̇ cos ϕ K r ot13 Δx cos ϕ cos η + K r ot22 sin ϕ
Iy
Conflict of interest The authors declare that they have no con-
flict of interest. 2
+ − |ϕ̇| cos2 η sin ϕ (K 21 x h + K 31 sin ϕ)
Iy
K v μ1 f mx h
− cos2 ϕ − Zq ,
Appendix A Iy Iy
π
where K r otmn = πρ( 21 − Δx̂)Imn and K v = 16 ρ I04 .
The higher-order dependence of the non-autonomous
The hinge line is set at 30% c (Δx̂ = 0.05), and the
aerodynamic vector field ga (x, t) (see Eqs. (59) and
value of C L α is computed based on the aspect ratio of
(60)) on the state vector x can be neglected. Thus, we
the wing, utilizing the extended lifting theory [70]. This
retain only the linear terms and obtain
flight dynamic model has been developed in Ref. [50],
⎛ ⎞ ⎛ ⎞ ⎛ 1 ⎞
u̇ −qw − g sin θ m X 0 (t) and the resulting eigenvalues of the averaged, linearized
⎜ ẇ ⎟ ⎜ qu + g cos θ ⎟ ⎜ ⎟
m Z 0 (t) ⎟
1
⎜ ⎟=⎜ ⎟+⎜ dynamics have been validated against numerical sim-
⎝ q̇ ⎠ ⎝ 0 ⎠ ⎝ M0 (t) ⎟
⎜ 1
⎠ ulations of Navier–Stokes equations by Sun et al. [55]
Iy
θ̇ q 0 and the experimental data of Cheng and Deng [59].
⎡ ⎤⎛ ⎞
X u (t) X w (t) X q (t) 0 u
⎢ Z u (t) Z w (t) Z q (t) 0 ⎥ ⎜ w ⎟
+⎢ ⎥⎜ ⎟
⎣ Mu (t) Mw (t) Mq (t) 0 ⎦ ⎝ q ⎠ . References
0 0 0 0 θ 1. Vela, P.A., Burdick, J.W.: Control of biomimetic locomo-
Assuming a horizontal stroke plane, parameterized by tion via averaging theory. In: IEEE International Conference
on Robotics and Automation, vol. 1, pp. 1482–1489. IEEE
the “back-and-forth” flapping angle ϕ, and a piecewise (2003)
constant variation in the wing pitch angle η, one obtains 2. Vela, P.A., Burdick, J.W.: Control of underactuated mechan-
[50]: X 0 (t) = −2K 21 ϕ̇(t)|ϕ̇(t)| cos ϕ(t) sin2 η, Z 0 (t) ical systems with drift using higher-order averaging theory.
= −K 21 ϕ̇(t)|ϕ̇(t)| sin 2η, M0 (t) = 2ϕ̇(t)|ϕ̇(t)| sin η In: Proceedings of the 42nd IEEE Conference on Decision
and Control, vol. 3, pp. 3111–3117. IEEE (2003)
[K 22 Δx̂ cos ϕ(t)+K 21 x h cos η(t)+K 31 sin ϕ(t) cos η], 3. Morgansen, K.A., Duidam, V., Mason, R.J., Burdick,
where x h is the distance from the vehicle center of J.W., Murray, R.M.: Nonlinear control methods for planar
mass to the root of the wing hinge line (i.e., the inter- carangiform robot fish locomotion. In: Proceedings 2001
section of the hinge line with the xb -axis) and Δx̂ ICRA. IEEE International Conference on Robotics and
Automation. vol. 1, pp. 427–434. IEEE (2001)
is the chord-wise distance from the center of pres- 4. Morgansen, K.A., Vela, P.A., Burdick, J.W.: Trajectory sta-
sure to this same hinge location, normalized by the bilization for a planar carangiform robot fish. In: Proceed-
chord length. Also, ρ is the air density, C L α is the ings. ICRA’02. IEEE International Conference on Robotics
three-dimensional lift curve slope of the wing, c(r ) is and Automation, vol. 1, pp. 756–762. IEEE (2002)
5. Wood, R.J.: The first takeoff of a biologically inspired at-
the spanwise
 R chord distribution, R is the wing radius, scale robotic insect. IEEE Trans. Robot. Autom. 24(2), 341–
Imn = 2 0 r m cn (r ) dr , and K mn = 41 ρC L α Imn . The 347 (2008)
time-varying stability derivatives are written directly 6. Taha, H., Hajj, M.R., Nayfeh, A.H.: Flight dynamics and
in terms of the system parameters as in [50]: X u = control of flapping-wing mavs: a review. Nonlinear Dyn.
70(2), 907–939 (2012)
−4 Km11 |ϕ̇| cos2 ϕ sin2 η, X w = − Km11 |ϕ̇| cos ϕ sin 2η, 7. Bullo, F.: Averaging and vibrational control of mechanical
X q = Km21 |ϕ̇| sin ϕ cos ϕ sin 2η − x h X w , Z u = 2X w , systems. SIAM J. Control Optim. 41(2), 542–562 (2002)
8. Taha, H.E., Woolsey, C.A., Hajj, M.R.: Geometric control
Z w = −2 Km11 |ϕ̇| cos2 η, Z q = 2 Km21 |ϕ̇| sin ϕ cos2 η −
approach to longitudinal stability of flapping flight. J. Guid.
K r ot12 K 12 Δx
m ϕ̇ cos ϕ − x h Z w , Mu = 4 I y |ϕ̇| cos ϕ sin η +
2 Control Dyn. 39(2), 214–226 (2015)
 K 12 Δx 9. Tahmasian, S., Woolsey, C.A.: A control design method
I y 2X q − x h Z u , Mw = 2 I y |ϕ̇| cos ϕ cos η +
m
for underactuated mechanical systems using high-frequency
2 KI21
y
|ϕ̇| sin ϕ cos2 η − mx
I y Z w, and
h inputs. J. Dyn. Syst. Meas. Control 137(7), 071,004 (2015)
10. Meerkov, S.: Principle of vibrational control: theory and
applications. IEEE Trans. Autom. Control 25(4), 755–762
(1980)

123
Reconciliation of two higher-order averaging techniques

11. Baillieul, J., Lehman, B.: Open-loop control using oscilla- 35. Nijmeijer, H., Van der Schaft, A.: Nonlinear Dynamical
tory inputs. In: CRC Control Handbook, pp. 967–980 (1996) Control Systems, vol. 175. Springer, Berlin (1990)
12. Markus, L., Yamabe, H.: Global stability criteria for differ- 36. Isidori, A.: Nonlinear Control Systems. Springer, Berlin
ential systems. Osaka Math. J. 12(2), 305–317 (1960) (2013)
13. Nayfeh, A.H., Mook, D.T.: Nonlinear Oscillations. Wiley, 37. Sastry, S.: Nonlinear Systems: Analysis, Stability, and Con-
Hoboken (2008) trol, vol. 10. Springer, Berlin (2013)
14. Nayfeh, A.H., Balachandran, B.: Applied Nonlinear 38. Sarychev, A.: Stability criteria for time-periodic systems via
Dynamics: Analytical, Computational, and Experimental high-order averaging techniques. In: Nonlinear Control in
Methods. Wiley, Hoboken (2008) the Year 2000, Lecture Notes in Control and Information
15. Nayfeh, A.H.: Perturbation Methods. Wiley, Hoboken Sciences, vol. 2, pp. 365–377. Springer (2001)
(2008) 39. Agrachev, A.A., Gamkrelidze, R., Sarychev, A.: Local
16. Khalil, H.K.: Nonlinear Systems. Prentice-Hall, New Jersey invariants of smooth control systems. Acta Appl. Math.
2(5), 1–5 (1996) 14(3), 191–237 (1989)
17. Sanders, J.A., Verhulst, F., Murdock, J.: Averaging: the peri- 40. Grobner, W., Knapp, H.: Contributions to the method of Lie
odic case. In: Averaging Methods in Nonlinear Dynamical series (1967)
Systems, pp. 21–44. Springer (2007) 41. Hale, J.: Ordinary Differential Equations. John Wiley, New
18. Murdock, J.A.: Perturbations: Theory and Methods, vol. 27. York (1969)
SIAM, Philadelphia (1999) 42. Mathematica file associated with the third and fourth order
19. Guckenheimer, J., Holmes, P.: Nonlinear Oscillations, proof. http://taha.eng.uci.edu/biblio-html/index.html
Dynamical Systems, and Bifurcations of Vector Fields, vol. 43. Kapitza, P.: Dynamic stability of a pendulum with an oscil-
42. Springer, Berlin (2013) lating point of suspension. J. Exp. Theor. Phys. 21(5), 588–
20. Berg, J.M., Wickramasinghe, I.M.: Vibrational control with- 597 (1951)
out averaging. Automatica 58, 72–81 (2015) 44. Taha, H.E., Tahmasian, S., Woolsey, C.A., Nayfeh, A.H.,
21. Bullo, F., Lewis, A.D.: Geometric control of mechanical M.R, Hajj: The need for higher-order averaging in the sta-
systems. In: Applied Mathematics. Springer, Berlin (2004) bility analysis of hovering, flapping-wing flight. Bioinspir.
22. Hassan, A.M., Taha, H.E.: Combined averaging-shooting Biomim. 10(1), 016,002 (2015)
approach for the analysis of flapping flight dynamics. J. 45. Stephenson, A.: Xx. On induced stability. Lond. Edinb.
Guid. Control Dyn. 41(2), 542–549 (2017) Dublin Philos. Mag. J. Sci. 15(86), 233–236 (1908)
23. Hassan, A.M., Taha, H.E.: Differential-geometric-control 46. Kapitza, P.L.: Dynamical stability of a pendulum when its
formulation for flapping flight multi-body dynamics. J. Non- point of suspension vibrates, and pendulum with a vibrating
linear Sci. 1–39 (2018) suspension. Collect. Pap. PL Kapitza 2, 714–737 (1965)
24. Sanders, J.A., Verhulst, F.: Averaging Methods in Nonlinear 47. Stabilization of an inverted pendulum under high-frequency
Dynamical Systems. Springer, Berlin (1985) excitation (Kapitza pendulum). https://www.youtube.com/
25. Bogoliubov, N.N., Mitropolsky, Y.A., Gillis, J.: Asymptotic watch?v=is_ejYsvAjY
methods in the theory of non-linear oscillations. Phys. Today 48. Sun, M.: Insect flight dynamics: stability and control. Rev.
16, 61 (1963) Mod. Phys. 86(2), 615 (2014)
26. Mitropolsky, I.A.: Averaging method in non-linear mechan- 49. Nelson, R.C.: Flight Stability and Automatic Control.
ics. Int. J. Non-linear Mech. 2(1), 69–96 (1967) McGraw-Hill, New York (1989)
27. Vela, P.A.: Averaging and Control of Nonlinear Systems. 50. Taha, H.E., Hajj, M.R., Nayfeh, A.H.: On the longitudinal
Ph.D. Thesis, California Institute of Technology (2003) flight dynamics of hovering MAVs/insects. J. Guid. Control
28. Agrachev, A., Gamkrelidze, R.: Chronological algebras and Dyn. 37(3), 970–978 (2014)
nonstationary vector fields. J. Math. Sci. 17(1), 1650–1675 51. Alexander, D.E.: Nature’s Flyers: Birds, Insects, and the
(1981) Biomechanics of Flight. JHU Press (2002)
29. Agračev, A., Gamkrelidze, R.V.: The exponential represen- 52. Taylor, G.K., Thomas, A.L.R.: Animal flight dynamics II.
tation of flows and the chronological calculus. Math. USSR- Longitudinal stability in flapping flight. J. Theor. Biol. 214,
Sbornik 35(6), 727 (1979) 351–370 (2002)
30. Fliess, M., Lamnabhi, M., Lamnabhi-Lagarrigue, F.: An 53. Taylor, G.K., Thomas, A.L.: Dynamic flight stability in the
algebraic approach to nonlinear functional expansions. IEEE desert locust schistocerca gregaria. J. Exp. Biol. 206(16),
Trans. Circuits Syst. 30(8), 554–570 (1983) 2803–2829 (2003)
31. Sarychev, A.: Stability criteria for time-periodic systems via 54. Sun, M., Xiong, Y.: Dynamic flight stability of a hovering
high-order averaging techniques. In: Nonlinear Control in bumblebee. J. Exp. Biol. 208(3), 447–459 (2005)
the Year 2000, vol. 2, pp. 365–377. Springer (2001) 55. Sun, M., Wang, J., Xiong, Y.: Dynamic flight stability of
32. Nayfeh, A.H.: Introduction to Perturbation Techniques. hovering insects. Acta Mech. Sin. 23(3), 231–246 (2007)
Wiley, Hoboken (1981) 56. Xiong, Y., Sun, M.: Dynamic flight stability of a bumblebee
33. Yagasaki, K., Ichikawa, T.: Higher-order averaging for peri- in forward flight. Acta Mech. Sin. 24(1), 25–36 (2008)
odically forced weakly nonlinear systems. Int. J. Bifurc. 57. Gao, N., Aono, H., Liu, H.: A numerical analysis of dynamic
Chaos 9(03), 519–531 (1999) flight stability of hawkmoth hovering. J. Biomech. Sci. Eng.
34. Sarychev, A.V.: Lie-and chronologico-algebraic tools for 4(1), 105–116 (2009)
studying stability of time-varying systems. Syst. Control 58. Faruque, I., Humbert, J.S.: Dipteran insect flight dynam-
Lett. 43(1), 59–76 (2001) ics. Part 1. Longitudinal motion about hover. J. Theor. Biol
264(2), 538–552 (2010)

123
 M. Maggia et al.

59. Cheng, B., Deng, X.: Translational and rotational damping 66. Popovic, P., Nayfeh, A.H., Oh, K., Nayfeh, S.A.: An
of flapping flight and its dynamics and stability at hovering. experimental investigation of energy transfer from a high-
IEEE Trans. Robot. 27(5), 849–864 (2011) frequency mode to a low-frequency mode in a flexible struc-
60. Khan, Z.A., Agrawal, S.K.: Control of longitudinal flight ture. Modal Anal. 1(1), 115–128 (1995)
dynamics of a flapping-wing micro air vehicle using time- 67. Oh, K., Nayfeh, A.: High-to low-frequency modal interac-
averaged model and differential flatness based controller. In: tions in a cantilever composite plate. J. Vib. Acoust. 120(2),
American Control Conference, 2007. ACC’07, pp. 5284– 579–587 (1998)
5289. IEEE (2007) 68. Taha, H.E., Nayfeh, A.H., Hajj, M.R.: Effect of the
61. Schenato, L., Campolo, D., Sastry, S.: Controllability issues aerodynamic-induced parametric excitation on the lon-
in flapping flight for biomimetic micro aerial vehicles gitudinal stability of hovering MAVs/insects. Nonlin-
(MAVs). In: Proceedings of the 42nd IEEE Conference on ear Dyn. 78, 2399–2408 (2014). https://doi.org/10.1007/
Decision and Control, vol. 6, pp. 6441–6447. IEEE (2003) s11071-014-1596-6
62. Nayfeh, S., Nayfeh, A.: Nonlinear interactions between two 69. Taha, H.E., Kiani, M., Hedrick, T.L., Greeter, J.S.M.: Vibra-
widely spaced modes—external excitation. Int. J. Bifurc. tional control: A hidden stabilization mechanism in insect
Chaos 3(02), 417–427 (1993) flight (under review). Nat. Commun
63. Nayfeh, S., Nayfeh, A.: Energy transfer from high-to low- 70. Taha, H.E., Hajj, M.R., Beran, P.S.: Unsteady nonlinear
frequency modes in a flexible structure via modulation. J. aerodynamics of hovering MAVs/insects. In: AIAA-Paper
Vib. Acoust. 116(2), 203–207 (1994) 2013–0504 (2013)
64. Nayfeh, A., Mook, D.: Energy transfer from high-frequency
to low-frequency modes in structures. J. Vib. Acoust.
117(B), 186–195 (1995) Publisher’s Note Springer Nature remains neutral with regard
65. Nayfeh, A.H., Chin, C.M.: Nonlinear interactions in a para- to jurisdictional claims in published maps and institutional affil-
metrically excited system with widely spaced frequencies. iations.
Nonlinear Dyn. 7(2), 195–216 (1995)

123