Identification of a Hybrid Spring Mass Damper via Harmonic Transfer Functions as a Step Towards Data-Driven Models for Legged Locomotion

0.6 MB·English

by  İsmail Uyanık

Identification of a Hybrid Spring Mass Damper via
Harmonic Transfer Functions as a Step Towards
Data-Driven Models for Legged Locomotion
˙Ismail Uyanık∗, M. Mert Ankaralı†, Noah J. Cowan†, O¨mer Morgu¨l∗ and Uluc¸ Saranlı‡
∗Dept. of Electrical and Electronics Engineering, Bilkent University, 06800 Ankara, Turkey
Email: [email protected], [email protected]
5
1 †Dept. of Mechanical Engineering, Johns Hopkins University, Baltimore, MD 21218, USA
0 Email: [email protected], [email protected]
2 ‡Dept. of Computer Engineering, Middle East Technical University, 06800 Ankara, Turkey

n Email: [email protected]
a
J Abstract—There are limitations on the extent to which man- construction of approximate solutions to such non-integrable

2 ually constructed mathematical models can capture relevant hybrid models [6–9].
2 aspects of legged locomotion. Even simple models for basic
behaviours such as running involve non-integrable dynamics, When accurate analytical solutions to the dynamics of
] requiring the use of possibly inaccurate approximations in the a legged platform are available [8], their structure can be
O design of model-based controllers. In this study, we show how exploited to yield effective solutions for system identifica-
R data-driven frequency domain system identification methods can tion and adaptive control [10]. Despite our previous studies
be used to obtain input–output characteristics for a class of dy- showing how accurate such models may be, there will always
.
s namicalsystemsaroundtheirlimitcycles,withhybridstructural be unmodeled components in the physical system, resulting
c propertiessimilartothoseobservedinleggedlocomotionsystems.
in discrepancies between the model and experiments [11].
[
Undercertainassumptions,wecanapproximatehybriddynamics
Attempts to manually incorporate these effects into the model
of such systems around their limit cycle as a piecewise smooth
1 is daunting at best, and often impossible. Consequently, we
linear time periodic system (LTP), further approximated as a
v proposeanalternativemethodinthisstudy,namelyusingdata-
time-periodic,piecewiseLTIsystemtoreduceparametricdegrees
8
of freedom in the identification process. In this paper, we use a drivensystemidentificationmethodstoderiveaninput–output
2
simple one-dimensional hybrid model in which a limit-cycle is transferfunctionforsuchhybridleggedlocomotionbehaviors,
6
induced through the actions of a linear actuator to illustrate thereby eliminating the need to manually construct an explicit
5
the details of our method. We first derive theoretical harmonic mathematical model for the system.
0
transfer functions of our example model. We then excite the
.
1 modelwithsmallchirpsignalstointroduceperturbationsaround Our main goal in this study is to provide a system
0 its limit-cycle and present systematic identification results to identification framework applicable to a useful (although not
5 estimate the harmonic transfer functions for this model. Com- comprehensive) class of legged locomotion models [8], and
1 parison between the data-driven HTF model and its theoretical possibly more complex robotic systems [12]. Our approach
: predictionillustratesthepotentialeffectivenessofsuchempirical is based on considering legged locomotion as a hybrid non-
v
i identification methods in legged locomotion. linear dynamical system with a stable periodic orbit (limit-
X
cycle), corresponding to the locomotor behavior of interest.
r I. INTRODUCTION We introduce a formulation that addresses the input–output
a
system identification problem in the frequency domain for a
Legged locomotion emerges from a staggering diversity
sub-class of hybrid legged locomotion models. More specifi-
of animal and robot morphologies and gaits, and modeling
cally, following certain assumptions on the hybrid dynamics
locomotordynamicsremainsagrandchallengeinbothbiology
of legged systems, we approximate their hybrid dynamics
and robotics [1,2]. Running behaviors, in particular, are com-
around the limit-cycle as a linear time-periodic system (LTP).
monly represented by relatively simple spring–mass models
However, this first LTP approximation is infinite dimensional,
such as the Spring-Loaded Inverted Pendulum (SLIP) model
makingparametricidentificationchallenging.Wehencefurther
[3,4]. A common feature of such models, however, is that
approximate the dynamics as a finite dimensional piecewise
their hybrid system dynamics involve intermittent foot contact
LTI system (maintaining its LTP nature), thereby limiting
with the ground, alternating between flight and stance phases
the parametric degrees of freedom while enabling a practical
during locomotion. Despite the presence of seemingly simple
identification framework.
models for basic behaviors such as running and walking, the
hybriddynamicsassociatedwiththesebehaviorscanberather Existing studies on system identification of LTP systems
complex,withnon-integrablepartssuchasthestancephase[3, focus on modeling these systems as multi-input single-output
5]. Given the utility of having accurate models and associated LTI systems. This approach is based on the concept of har-
analyticsolutionsinconstructinghighperformancecontrollers monic transfer functions [13], which are infinite-dimensional
fornonlinearsystems,substantialefforthasbeendevotedtothe operators that are analogous to frequency response functions
for LTI systems. An identification strategy for such systems q¯(t) = q¯(t−T) when u(t) = 0. For such systems, if we let
was developed in [14] using power spectral density and cross q(t) = q¯(t)+x(t) and linearize the dynamics in (1) around
spectral density functions. A similar method was used in [15] the limit-cycle q¯(t), and u(t)=0 we get
considering the effects of noise in both input and output
x˙(t)=A(t)x(t)+B(t)u(t)
measurements. Different than these studies, local polynomial (2)
methods and lifting approaches were also used for the identi- y(t)=C(t)x(t)+D(t)u(t)
fication of harmonic transfer functions for multi-input single- where
output models of LTP systems [16]. (cid:20) (cid:21)
∂f
A(t) = , (3)
Our contributions in this paper focus on representing the ∂q q(t) =q¯(t)
dynamics of legged locomotion as a linear time periodic u(t)= 0
system, thereby enabling the use of the system identification (cid:20)∂f(cid:21)
B(t) = . (4)
method proposed in [14] for such systems. We achieve this ∂u q(t) =q¯(t)
by using a new phase definition in identifying the harmonic u(t)= 0
transfer functions, illustrated in the context of a simplified
ThiscorrespondstoaLinearTimePeriodic(LTP)system,with
model designed to mirror structural properties of legged lo-
all system matrices sharing a common period, T.
comotion models. When the problem is approached as a grey-
box model with finite parameters (piecewise LTI), it suffices
B. Modeling Framework for Hybrid Systems
to non-parametrically estimate a finite number of harmonics,
to which we later fit parametric models. Legged systems are often modeled using hybrid dynamics
due to intermittent foot contact with the ground, which cannot
Ankarali and Cowan [17] developed a similar system
be represented with a single, smooth dynamical flow. In
identification method for hybrid systems with periodic orbits
the broadest sense, a hybrid dynamical system is a set of
using “discrete time” harmonic transfer functions. However,
smoothflowstogetherwithdiscretetransitions(andassociated
the framework and assumptions in this paper are distinctly
transformations)betweentheseflowstriggeredbyintersections
different from their approach. Specifically, they use mappings
of system trajectories with sub-manifolds of the continuous
between different cross sections to construct a discrete-time
state space [21]. These flows are called charts, indexed with
LTP system, and use discrete time HTFs for identification.
unique labels I := {0,··· ,d} each with possibly different
Also, the current paper focuses on harmonic balance, which
equations of motion. Along its trajectories, a hybrid system
also distinguishes this paper from [17].
transitions from one chart to another, with each transition
defined by the zero crossing of a threshold function. For each
II. REPRESENTATIONOFLEGGEDLOCOMOTIONASA source chart α∈I and destination chart β ∈I, the threshold
HYBRIDDYNAMICALSYSTEM function hβ defines the transition from chart α to chart β.
α
Our goal in this study is to provide a system identification An example transition graph for a hybrid dynamical system is
framework for a class of models related to legged locomotion illustrated in Fig. 1.
using harmonic transfer functions. For the present paper,
Since we are interested in the local behaviour around
we limit ourselves to “clock-driven” locomotion models as
the limit-cycle, we assume that there is only one transition
described in Section II-A, representative of controllers used function associated with each chart.1 We further assume that
by a wide variety of existing robots [12,18,19], with open-
system trajectories are continuous at transitions, meaning that
loop central pattern generators (CPG) coordinating control
system states do not experience discrete changes coincident
actions to achieve time periodic behaviour. This will allow
with chart transitions. As a final note, we assume that the
us to directly use time periodicity in our LTP analysis, while
hybrid dynamical system we consider has an isolated periodic
eliminatingavarietyofcomplicationsassociatedwithestimat-
orbit ensuring that chart transitions within the limit cycle are
ing the phase [20].
also periodic and consistent.
A. Smooth Clock-driven Oscillators It is important to note that these assumptions are generally
satisfied by models of common locomotory behaviors such as
Ingeneral,thedynamicsofsmooth,clock-drivenoscillators running and walking [8,22] as well as a wide range of legged
with external inputs can be written as robots for which leg masses are negligible compared to the
q˙=f(q,φ,u),
1This approach does not apply to gaits such as pronking that nominally
φ˙ =1 involvemultiplelegsmakingcontactwiththegroundatthesametimewhen
(1) onthelimitcycle,becausesmalldeviationsfromthelimitcyclecanleadto
f :Rn×S1×Rp (cid:55)→Rn
differenttouch-downorderbetweenlegs,violatingourassumption.
(q, φ)∈ Rn×S1, u∈Rq
where (q,φ) and Rn × S1 denote the state vector and the hβα
statespaceoftheoscillator,respectively.Thecirclecomponent
S1 = mod(R+,T) enforces the periodicity of the dynamics, Chartα Chartβ
while the external input u(t) represents small external pertur-
bations which we will use for system identification. hα
β
In this paper, we focus on oscillators of the form (1)
Fig.1. Asimplestatetransitiongraphforahybriddynamicalsystem.
withasymptoticallystable,isolatedperiodicorbits(limit-cycle)
dynamics of a larger body [12,18]. Consequently, the system the LTP structure of the system. The LTP equations of motion
identification methods we introduce will remain applicable to then take the form
systems other than the simplified example we will present in (cid:40)
A x(t)+B u(t), if mod(t,T)∈[0,tˆ)
this paper. x˙(t)≈ 0 0 (8)
A x(t)+B u(t), if mod(t,T)∈[tˆ,T)
1 1
(cid:40)
C. Modeling Legged Locomotion as a Linear Time Periodic C0x(t)+D0u(t), if mod(t,T)∈[0,tˆ)
y(t)≈ (9)
System C x(t)+D u(t), if mod(t,T)∈[tˆ,T)
1 1
Forclarity,welimitourfocusinthissectiontoanexample The formulation above constitutes the basis of our framework
hybrid dynamical system with only two charts, I = {0,1}, for analyzing and identifying clock-driven legged locomotion
designed to capture stance and flight phases of simple spring- models.
mass models of locomotion. Based on a clock driven as-
sumption, for each i ∈ I the continuous dynamics can be
III. HARMONICTRANSFERFUNCTIONS
represented with
A. Preliminaries and Background
φ˙ =1
System identification studies on (stable) LTI systems rely
q˙ =f (q,φ,u) , (5)
i i on the fact that if the input is sinusoidal, then, at steady-
qi ∈Rn state, the output will also be a sinusoidal signal (at the
same frequency but with a possibly different magnitude and
and let the associated threshold function be hmod(i+1,2)(q). phase). This one-to-one mapping between input and output
i
Thetransitionmapassociatedwitheachhybrideventissimply signals allows us to characterize the dynamics in terms of a
the identity map, q (cid:55)→ q , due to the continuity assumption. frequency response function (FRF) also know as a Bode plot.
i i
Our linearization of these hybrid dynamics towards an LTP Unfortunately, this approach does not readily transfer to LTP
approximation assumes that these transition times, tˆ, zero systems, which produce output spectra that include multiple
crossings of h1(q) and h0(q), maintain their periodicity and (possibly infinite) harmonics of the input stimuli, each with
0 1
offsets within the period in close proximity of the limit-cycle, possibly different magnitude and phase at steady state.
resulting in the following form of the nonlinear dynamics
One ad hoc way to mitigate this is to enforce a one-to-one
φ˙ =1 (6) mapping by neglecting higher harmonics [23]. However, this
(cid:40) assumption may result in substantial inaccuracies particularly
f0(q,φ,u) , if mod(t,T)∈[0,tˆ) when the influence of higher harmonics on the response is
q˙≈ . (7)
f (q,φ,u) , if mod(t,T)∈[tˆ,T) expectedtobesignificant.Motivatedbythisproblem,Wereley
1
[24] proposed a linear one-to-one mapping between the coef-
Assuming that the system given above has a limit cycle q¯(t) ficients of an exponentially modulated periodic (EMP) signal
withaperiodT,linearizationaroundq¯(t)yieldsthepiecewise at the input of LTP systems to the coefficients of an EMP
smooth LTP system signal at their output. This linear operator that maps the input
harmonics to the output harmonics of an LTP system is called
(cid:40)
A (t)x(t)+B (t)u(t), if mod(t,T)∈[0,tˆ) a Harmonic Transfer Function (HTF) [13].
0 0
x˙(t)=
A1(t)x(t)+B1(t)u(t), if mod(t,T)∈[tˆ,T) In the following section, we review the derivation of
harmonictransferfunctionsaspresentedin[13]and[25],using
where theprincipleofharmonicbalancestartingfromthestatespace
(cid:20)∂f (cid:21) (cid:20)∂f (cid:21) representation of (2).
A (t):= 0 ,B (t):= 0 ,
0 ∂q q(t) =q¯(t) 0 ∂u q(t) =q¯(t)
u(t)= 0 u(t)= 0 B. Theoretical Derivation of Harmonic Transfer Functions
(cid:20) (cid:21) (cid:20) (cid:21)
∂f ∂f
A (t):= 1 ,B (t):= 1 . Recall that the system matrices in (2) are all T-periodic.
1 ∂q q(t) =q¯(t) 1 ∂u q(t) =q¯(t) Consequently, they can be represented by an infinite Fourier
u(t)= 0 u(t)= 0 series with pumping frequency ω = 2π/T. For the system
p
matrix A(t), we have
It is natural to assume that direct measurement of all x(t)
∞
maynotbeavailableorwemayonlymeasureasubsetofx(t). (cid:88)
A(t)= A ejωpnt . (10)
Consequently, we also define a time-periodic output equation n
n=−∞
as in the form (9).
The matrices B(t), C(t) and D(t) can be similarly decom-
Since system matrices A (t), B (t), C (t) and D (t) with
i i i i posed. In addition, we can also expand the state and output
i ∈ {0,1} are time parametrized functions, the system has
vectors since they are both EMP signals. Substituting these
infiniteparametricdegreesoffreedom,makingparametricsys-
expansions into (2) and applying the principle of harmonic
tem identification challenging even when Harmonic Transfer
balance as explained in [13], we obtain the harmonic state
Functions are used. At this point, we hypothesize that for
space representation as
hybridsystems,thevariabilitywithinachartissmallcompared
to the change between charts and we approximate the LTP sX =(A−N)X +BU
(11)
dynamics using a piecewise LTI approximation that preserves Y =CX +DU,
where the doubly infinite vectors representing the harmonics C. Estimation of Harmonic Transfer Functions via Frequency
of the state, control, and output signals are Domain System Identification
XT :=[··· ,xT ,xT ,xT,xT,xT,···], In this section, we briefly explain the data-driven system
−2 −1 0 1 2
identification method presented in [14].
UT :=[··· ,uT ,uT ,uT,uT,uT,···], (12)
−2 −1 0 1 2
InanLTPsystem,asinusoidalinputataspecificfrequency
YT :=[··· ,yT ,yT ,yT,yT,yT,···],
−2 −1 0 1 2 generates a superposition of sinusoids at multiple (possibly
and the doubly infinite input modulation matrix is an infinite number of) harmonics. Consequently, the system
identification framework starts with truncating the number of
N :=blockdiag{jnωpI}, ∀n∈Z, (13) harmonictransferfunctionsGˆ tobeestimated.Inthefollowing
examples, we consider only three frequencies in the output to
which modulates the input frequency to different harmonic
clearly illustrate the approach of [14].
frequencies. Details on the derivations can be found in [13].
Suppose that an LTP system consists of the superposition
The T-periodic dynamics matrix, A(t), is expressed in of three different harmonic transfer functions, Gˆ , Gˆ and
terms of its complex Fourier coefficients, {An|n ∈ Z}, as Gˆ ,eachcorrespondingtoadifferentfrequencyco0mpo−ne1ntof
1
a doubly infinite block Toeplitz matrix,
the output. The output can then be expressed as
... ... ... ... ... ...  Yˆ(jω) = Gˆ0(jω)U(jω)+Gˆ−1(jω)U(jω+jωp)
··· A0 A−1 A−2 A−3 A−4 ··· +Gˆ1(jω)U(jω−jωp). (17)
··· A A A A A ···
 1 0 −1 −2 −3 
A=··· A A A A A ···, (14) Inthisnewformulation,thenth transferfunctionisdefined
 2 1 0 −1 −2 
··· A A A A A ··· asthelinearoperatorthatmapstheoutputatfrequencyωtoan
 3 2 1 0 −1 
··· A4 A3 A2 A1 A0 ··· input at the same frequency,modulated with ejnωpt. However,
 
... ... ... ... ... ... a single input–output pair in each frequency will naturally
not be sufficient to estimate harmonic transfer functions as
in the case of LTI systems, since the identification problem
with a similar definition for B(t) in terms of its Fourier
coefficients represented by {B |n ∈ Z}, C(t) in terms of will then be underdetermined. Therefore, either at least three
{C |n∈Z}, and D(t) in termsnof {D |n∈Z}. independent inputs or additional constraints must be provided
n n toenableasuccessfulidentificationoftheseharmonictransfer
This collection of doubly infinite matrices is called the functions.
harmonicstatespacemodel(HSS)ofthesystemgivenin(2).
There are two key issues that need to be addressed before
However, it will also be useful to determine an explicit input–
designing input signals for the identification process. First,
output functional relationship between the Fourier coefficients
we will require the use of at least as many variations on the
of the harmonics of the input, {u |n ∈ Z}, and those of
n input signal as the number of harmonic transfer functions to
the output, {y |n ∈ Z}. This relationship is represented by
n beestimated.Thisisaccomplishedin[14],whichusesasingle
the harmonic transfer function, G(s), which is also an infinite
input sequence signal for systemidentification, constructed by
dimensional matrix of Fourier coefficients, satisfying
concatenatingphaseshiftedcopiesofasinglewaveformonthe
input evenly separated by delays within the system period. A
Y =GU . (15)
completecharacterizationofsystemdynamicsispossiblewith
Based on (11), G can be computed as thismethodsincedifferentmodesofthesystemwereactivated
through the use of phase-shifted copies of a single waveform.
G=C[sI−(A−N)]−1B+D , (16)
Thesecondissueistheneedtoexciteallfrequencycompo-
as long as the inverse within this equation exists. nents within the system by providing input signals with a suf-
ficiently wide frequency spectrum. This can be accomplished
There are, however, two problems associated with the through the use of chirp signals, whose frequency varies with
harmonic transfer function as stated above. First, it is not time.Theuseofchirpinputsignals,combinedwiththeideaof
clear whether the inverse of the doubly infinite matrix in supplyingmultiple,phase-shiftedinputsequencesallowsusto
the definition of the harmonic transfer function will always obtainsufficientlyrichinput–outputdatatosupportthesystem
exist. This problem will be dealt with by an application of identification process.
the Floquet Theorem. Second, the harmonic transfer function
is a doubly infinite matrix operator, which cannot practically Using this data with input–output pairs, one can estimate
be implemented on a computer. This second problem will be theharmonictransferfunctionsofthesystem,sothattheerror
mitigatedbytruncatingtheHTFinordertoimplementanalysis betweenactualandestimatedoutputsisminimized.Therefore,
on a computer. we can convert the identification problem to an optimization
formulated as
Note that the theoretical definition of harmonic transfer
functions in [13], reviewed in this section, requires the state minimize J=(Y−UTGˆ)2. (18)
Gˆ
space representation of the system to be available. Our goal is
to estimate this theoretically computed transfer function G(s) However,notethatSiddiqi[14]combinesallphase-shifted
by using input-output data in the frequency domain without signals in a single input. Hence, the problem is still underde-
necessitating knowledge of internal system dynamics. termined in the frequency domain, since a single input–output
pair for a specific frequency will not be sufficient to identify Fig. 2 illustrates the vertical leg model we focus on in this
three harmonic transfer functions. In order to address this section. It consists of a mass attached to a leg with a spring-
problem, they consider additional constrains on the estimated damper mechanism as well as a force transducer. Unlike the
harmonic transfer functions. First, they assume that transfer SLIP model, we assume that the toe is permanently affixed
functionsaresmooth,whichisreasonableforphysicalsystems. on the ground. Nevertheless, we recover the hybrid nature of
This is enforced through a difference operator, D2, designed locomotory gaits by assuming that the damper is turned on
to compute the second derivative of a vector when multiplied duringa“stancephase”(lossy)andoffduringa“flightphase”
from the left side. Details on the derivations for D2 can be (lossless). This construction recovers the hybrid nature of the
found in [14]. dynamics,whileallowingactiveinputthroughouttheentiretra-
jectory to support the generation of system identification data,
The smoothness condition on transfer functions requires
as well as admitting theoretical computation of its harmonic
penalizingthecurvatureofindividualtransferfunctions.There-
transfer functions for a comparative investigation.
fore, [14] modifies (18) to include a cost associated with the
curvature,yieldingarevisedminimizationproblemformulated We use the force transducer f in this system for two
as purposes. Firstly, active energy input to the system must be
minimize J=(Y−UTGˆ)2+(αD2Gˆ)2 , (19) providedtomaintainthelimitcycleandcompensateforenergy
Gˆ losses due to the presence of damping. Second, it will be
where α is a manually tuned constant weight for penalizing used as an exogenous input to the system to support the
curvature. The solution of (19) can easily be obtained by system identification process. Many physical legged platforms
differentiating J with respect to Gˆ, taking the form include similar active components in their legs to regulate
their mechanical energy [26,27]. Notwithstanding differences
Gˆ =(UTU+αD4)−1UTY , (20) in how these actuators are incorporated into the system, they
can all be used as the necessary exogenous inputs to perform
where the rows of the matrix Gˆ(ω) correspond to individual
system identification. A similar model was also investigated
different harmonic transfer functions as
in [16] but using an additional nonlinear spring for energy
 Gˆ (ω)  regulation.
1
Gˆ(ω)= Gˆ0(ω) . (21) The equations of motion for this simplified legged loco-
Gˆ−1(ω) motion model are given by
(cid:26)
Note that UTU and UTY correspond to power spectral mx¨= −mg−cx˙ −k(x−x0)+f(t), if x˙ >0 (22)
and cross spectral density functions, respectively. Therefore, −g−k(x−x0)+f(t), otherwise.
(20) is analogous to estimating transfer functions in LTI
Thelossyandlosslessdynamicsin(22)correspondtodifferent
systems, with an additional cost on curvature.
charts in Fig. 1 and zero crossings of x˙ represent threshold
functions for both phases.
IV. SIMPLIFIEDLEGGEDLOCOMOTIONMODELWITH
HYBRIDSYSTEMDYNAMICS Ourillustrativeexamplesusetheparametersg =9.81,k =
200, c=2, m=1 and x =0.2, chosen to be similar to the
0
Inthissection,wedescribeasimple,verticallyconstrained
parametersofaverticalhopperplatforminourlaboratory[28].
spring-mass-damper system that possesses hybrid structural
As noted above, we choose the linear actuator input f(t) =
properties similar to the extensively studied Spring-Loaded
f (t)+u(t), consisting of a forcing term f (t) to compensate
0 0
Inverted Pendulum (SLIP) model for running behaviors. This
for energy losses, and a chirp signal u(t) to introduce small
will provide a simple example to illustrate the application of
periodic perturbations for system identification.
our system identification method to such systems.
B. Theoretical Computation of Harmonic Transfer Functions
A. System Dynamics
The goal of this section is to compute harmonic transfer
g functions for our model around its limit cycle as outlined in
M
Section III-B.
Wefirstassumethattheforcinginputf (t)isappropriately
0
chosen to induce an asymptotically stable limit cycle for this
system. For example, our simple leg model achieves a stable
k c f
limit cycle with f (t) = cos(2πt). At this point, changing
0
into error coordinates away from the limit cycle with ξ =
x(t)−x¯(t),andsubstitutinginto(22),theequationsofmotion
take the form
(cid:26)−cξ˙−kξ, if ξ˙+x¯˙(t)>0
ξ¨= (23)
−kξ, otherwise
Due to the simplicity of the dynamics, this corresponds to
Fig.2. Simplifiedlegmodel.
a piecewise LTI system without necessitating any additional
approximations, taking the form −20
(cid:20)ξ˙ (cid:21) (cid:20) 0 1 (cid:21)(cid:20)ξ (cid:21) (cid:20)0(cid:21) dB)
ξ˙21 = −k −cs(ξ˙,t) ξ12 + 1 u(t), (24) ude ( −40
nit
where the hybrid nature of the system is captured by the flag g
a
s(ξ˙,t), with s=1, when ξ˙+x¯˙(t)>0 and s=0 otherwise. M −60
We now need to represent this piecewise LTI system as a
linear time periodic system. However, even though the binary 0 Theoretica l
vliamluitedcyfuclnectiitosnelsf,(ξ˙t,hti)sciasnnboetctohnesicdaesreedfotirmter-apjeecritoodriiecsoanwtahye egree)−50 EPasrtaimmaetterdic
D
from the limit cycle. To proceed, we hence assume that input e (−100
induced perturbations are small, and that the binary valued as
function s(ξ˙,t) maintains its period and becomes strictly Ph−150
time dependent rather than state dependent, taking the form
s(ξ˙,t)≈s(t).WenowcanperformaFourierseriesexpansion 10 100
Frequency (Hz)
ons(t)bytreatingitasasquarewavewithanoffsettoobtain
a linear time periodic system in the form Fig.3. Estimationresultsforthefundamentalharmonic.
(cid:20)ξ˙ (cid:21) (cid:20) 0 1 (cid:21)(cid:20)ξ (cid:21) (cid:20)0(cid:21)
1 = 1 + u(t), (25)
ξ˙2 −k −cs(t) ξ2 1 Fig. 3 illustrates the estimation performance of our al-
(cid:20) (cid:21)
y = [1 0] ξ1 . gorithm for the magnitude and phase of the fundamental
ξ2 harmonic.Bothgraphsshowthattheapplicationoftheidentifi-
cationalgorithmin[14]workswellevenfornonlinearperiodic
PluggingtheseequationsintotheHTFframeworkdescribedin
systems with hybrid dynamics.
SectionIII-B,yieldsanalyticsolutionstotheharmonictransfer
functions. We omit the details of this derivation due to space Wealsoshowouridentificationresultsforthreeharmonics
considerations, but use the resulting analytic solutions for the in both the negative and positive sides in Fig. 4. Even though
harmonic transfer functions up to n = 10 to evaluate the magnitudes for the harmonic transfer functions are small
h
output of our system identification method. compared to the fundamental, the identification algorithm can
provide accurate estimates for these transfer functions except
C. Data-DrivenIdentificationofHarmonicTransferFunctions in some narrow regions of G−2 and G2. The identification
algorithm could not correctly estimate these two harmonics
In this section, we obtain harmonic transfer functions around12−15(rad/s).Onepossiblereasonforthisdiscrepancy
corresponding to the linearized dynamics of (25) by using is the presence of strong responses in all harmonics around
input–output data without assuming prior knowledge of the the same frequency except G and G , resulting in the
−2 2
state space model. Using f0(t) = cos(2πt) and u(t) = 0 for inability of the identification algorithm to distinguish between
30cycleswithoutaperturbation,ourexamplesystemstabilized thecontributionsfromeachharmonicabsentknowledgeofthe
to a limit cycle x¯(t) with a period T = 1s. We use the 30th internal system dynamics. Alternatively, these discrepancies
periodasthenumericallimitcycleofthenonlinearsystemand may also be a result of the fact that hybrid transitions are not
subtract it from the trajectories of subsequent experiments to strictly time periodic (rather, they are state-dependent) which
obtain the error function ξ1. likely has effects on different frequencies and harmonics. We
plan on investigating these issues further in the future.
In order to obtain input–output data for system identifica-
tion, we apply an input signal consisting of nine subsequent For a comparative analysis, we also present results from
30s long chirp signals, each with a linearly increasing fre- a parametric identification in order to show that further cor-
quency in the range (0,7] Hz over its duration but with a rections on estimation results from a non-parametric method
different starting phase evenly distributed across the system’s are possible. To this end, we fit the system parameters k and
period, T =1s. Each chirp signal has an amplitude of 0.004, c in (25) by comparing root mean square error between theo-
chosen to be large enough to perturb system dynamics but retically computed and estimated harmonic transfer functions
small enough to keep the system close to the periodic orbit. G , G and G . We truncate the system response after the
0 −1 1
A sample chirp signal with zero phase can be generated by first harmonic in order to discard erroneous regions in higher
u(t)=0.004sin(14πt2/30). (26) harmonics. The resulting estimates were kˆ = 200 for the
spring constant and cˆ = 2.12 for the damping coefficient,
which closely coincide with the parameters used to generate
The resulting output is then subtracted from the numeri-
the input–output data. As such, harmonic transfer functions
cally measured limit cycle to obtain error trajectories ξ for
1 obtained from parametric identification were found to closely
vertical position. The input signal and ξ are then used as
1 matchthoseobtainedfromtheoreticalcomputationsasseenin
in Section III-C to estimate harmonic transfer functions for
Fig. 4.
our system. Since our theoretical computations showed that
responsesbeyondthethirdharmonicwereverysmall,weonly Motivatedbytheseidentificationresults,weplantoextend
consider the fundamental harmonic and three harmonics on our work to the identification of the Spring-Loaded Inverted
both sides for our experiments. Pendulum(SLIP)model[3]anditsextensions,widelyusedas
xx 1100--33 x 10-3
3.2 5
gnitude 1.6 G-1 TEPhasertaiommreaettteircdiacl 2.5 G1
a
M
0 0
xx 1100--44 x 10-4
3.8 3
e
d G G
u -2 2
nit 1.9 1.5
g
a
M
0 0
xx 1100--33 x 10-3
1.4 7
ude G-3 G3
nit 0.7 3.5
g
a
M
0 0

5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
Frequency (rad/s) Frequency (rad/s)
Fig.4. Estimationresultsforthehigherorderharmonics.
models of locomotory behaviors in the literature. Our future functions of nonlinear periodic models with hybrid system
goalistoapplyoursystemidentificationmethodstoourphys- dynamics.
icalmonopodrobotplatformandtocomparetheidentification
performances with our previously verified analytical model
[11]. ACKNOWLEDGMENT
This material is based on work supported by the National
V. CONCLUSION Science Foundation (NSF) Grants 0845749 and 1230493 (to
N. J. Cowan). The authors thank Aselsan and The Scientific
In this paper, we presented a system identification strat-
and Technological Research Council of Turkey (TU¨B˙ITAK)
egy to estimate input–output transfer functions for a simple
for ˙Ismail Uyanık’s financial support.
hybrid spring mass damper system as a step towards data-
driven models for legged locomotion. We first showed that
a class of hybrid locomotion models can be approximated
REFERENCES
with a piecewise constant LTP systems in close proximity
to their asymptotically stable limit-cycle. Our analysis and [1] P.J.Holmes,R.J.Full,D.E.Koditschek,andJ.Guckenheimer,“The
dynamics of legged locomotion: Models, analyses, and challenges,”
identification framework is based on the concept of harmonic
SIAMRev,vol.48,no.2,pp.207–304,2006.
transfer functions [24].
[2] R.J.FullandD.E.Koditschek,“Templatesandanchors:neuromechan-
We first observed that the hybrid system dynamics associ- ical hypotheses of legged locomotion on land,” J Exp Biol, vol. 202,
no.23,pp.3325–3332,1999.
ated with this model exhibits piecewise LTI behavior around
its periodic orbit. We then represented this behavior as a [3] W. J. Schwind, “Spring loaded inverted pendulum running: A plant
model,” PhD Thesis, University of Michigan, Ann Arbor, MI, USA,
purely time periodic system around the limit cycle in order
1998.
to utilize system identification techniques applicable to Linear
[4] R. J. Full and M. S. Tu, “Mechanics of a rapid running insect: two-,
Time Periodic systems.
four-, and six-legged locomotion,” J Exp Biol, vol. 156, pp. 215–231,
1991.
In order to provide a basis for comparison, we computed
[5] P. Holmes, “Poincare´, celestial mechanics, dynamical-systems theory
analyticexpressionsforharmonictransferfunctionsassociated
and“chaos”.”PhysicsReports,vol.193,pp.137–163,September1990.
withtheLTPapproximationtooursimplifiedhybridmodel.In
[6] W.J.SchwindandD.E.Koditschek,“Approximatingthestancemap
our theoretical analysis, we considered the system’s response
of a 2 dof monoped runner,” Journal of Nonlinear Science, vol. 10,
up to the 10th harmonic. We observed that there were no no.5,pp.533–588,July2000.
meaningfulresponsesonbothpositiveandnegativesidesafter
[7] H.Geyer,A.Seyfarth,andR.Blickhan,“Spring-massrunning:Simple
the third harmonic. Therefore, we decided to truncate the sys- approximate solution and application to gait stability,” Journal of
temresponseafterthethirdharmonicduringouridentification TheoreticalBiology,vol.232,no.3,pp.315–328,February2005.
studies. [8] U.Saranli,O.Arslan,M.M.Ankaralı,andO¨.Morgu¨l,“Approximate
analytic solutions to non-symmetric stance trajectories of the passive
Wethenperformedsystematicsimulationstudiesandiden- spring-loadedinvertedpendulumwithdamping,”NonlinearDynamics,
tified the harmonic transfer functions of the same model vol.62,pp.729–742,December2010.
without knowledge of its internal dynamics. We used an [9] M. M. Ankarali and U. Saranli, “Stride-to-stride energy regulation
inputsignalconsistingofsuccessivechirpsignals,withphases for robust self-stability of a torque-actuated dissipative spring-mass
hopper,” Chaos: An Interdisciplinary Journal of Nonlinear Science,
evenly distributed across the system’s period, to obtain a full
vol.20,no.3,p.033121,September2010.
characterization of system dynamics for our frequency range
[10] I.Uyanik,U.Saranli,andO¨.Morgu¨l,“Adaptivecontrolofaspring-mass
of interest. Our studies showed that LTP system identification
hopper,”inIEEEInternationalConferenceonRoboticsandAutomation
techniques can successfully be used to identify the transfer (ICRA),2011,pp.2138–2143.
[11] I. Uyanik, O¨. Morgu¨l, and U. Saranli, “Experimental validation of a [20] S.RevzenandJ.M.Guckenheimer,“Estimatingthephaseofsynchro-
feed-forward predictor for the spring-loaded inverted pendulum tem- nizedoscillators,”PhysicalReview,vol.78,no.5,p.051907,November
plate,”IEEETransactionsonRobotics,2015. 2008.
[12] U. Saranli, M. Buehler, and D. E. Koditschek, “RHex: A simple [21] J. Guckenheimer and S. Johnson, “Planar hybrid systems,” in Hybrid
andhighlymobilerobot,”InternationalJournalofRoboticsResearch, systemsII. Springer,1995,pp.202–225.
vol.20,no.7,pp.616–631,July2001. [22] A.WuandH.Geyer,“The3-dspring–massmodelrevealsatime-based
[13] N. W. Wereley, “Analysis and control of linear periodically time deadbeat control for highly robust running and steering in uncertain
varying systems,” Ph.D., Massachusetts Institute of Technology, Dept. environments,”IEEETransactionsonRobotics,2013.
ofAeronauticsandAstronautics,1991. [23] A.Leonhard,“Thedescribingfunctionmethodappliedfortheinvesti-
[14] A. Siddiqi, “Identification of the harmonic transfer functions of a gationofparametricoscillations,”inProceedingsofthe2ndIFACWorld
helicopter rotor,” M.Sc., Massachusetts Institute of Technology, Dept. Conference,1963,pp.21–28.
ofAeronauticsandAstronautics,2001. [24] N. Wereley and S. Hall, “Frequency response of linear time periodic
[15] S.Hwang,“Frequencydomainsystemidentificationofhelicopterrotor systems,”inProceedingsofthe29thIEEEConferenceonDecisionand
dynamicsincorporatingmodelswithtimeperiodiccoefficients,”Ph.D., Control,1990,pp.3650–3655vol.6.
GraduateSchooloftheUniversityofMarylandatCollegePark,1997. [25] E.Mo¨llerstedt,“Dynamicanalysisofharmonicsinelectricalsystems,”
[16] M.W.SracicandM.S.Allen,“Methodforidentifyingmodelsofnon- Ph.D.,LundInstituteofTechnology,DepartmentofAutomaticControl,
linearsystemsusinglineartimeperiodicapproximations,”Mechanical 2000.
SystemsandSignalProcessing,vol.25,no.7,pp.2705–2721,2011. [26] G.PiovanandK.Byl,“Two-elementcontrolfortheactiveslipmodel,”
[17] M.M.AnkaraliandN.J.Cowan,“Systemidentificationofrhythmichy- inIEEEInternationalConferenceonRoboticsandAutomation(ICRA),
briddynamicalsystemsviadiscretetimeharmonictransferfunctions,” May2013,pp.5656–5662.
in Proc IEEE Int Conf on Decision and Control, Los Angeles, CA, [27] G.SecerandU.Saranli,“Controlofmonopedalrunningthroughtunable
USA,December2014. damping,” in Signal Processing and Communications Applications
[18] K.C.Galloway,G.C.Haynes,B.D.Ilhan,A.M.Johnson,R.Knopf, Conference(SIU),201321st,April2013,pp.1–4.
G. Lynch, B. Plotnick, M. White, and D. E. Koditschek, “X-rhex: A [28] I. Uyanik, “Adaptive control of a one-legged hopping robot through
highly mobile hexapedal robot for sensorimotor tasks,” University of dynamically embedded spring-loaded inverted pendulum template,”
Pennsylvania,Tech.Rep.,2010. M.Sc.,BilkentUniv.,Ankara,Turkey,August2011.
[19] U. Saranli, A. Avci, and M. C. O¨ztu¨rk, “A modular real-time field-
bus architecture for mobile robotic platforms,” IEEE Transactions on
InstrumentationandMeasurement,vol.60,no.3,pp.916–927,2011.