NASA/TM–2014–218314

Exploratory Study for

Continuous-Time Parameter

Estimation of Ankle Dynamics

Sunil L. Kukreja

NASA Armstrong Flight Research Center

Richard D. Boyle

NASA Ames Research Center

May 2014
NASA STI Program . . . in Profile

Since its founding, NASA has been

dedicated to the advancement of

aeronautics and space science. The NASA

scientific and technical information (STI)

program plays a key part in helping NASA

maintain this important role.

The NASA STI Program operates under

the auspices of the Agency Chief

Information Officer. It collects, organizes,

provides for archiving, and disseminates

NASA’s STI. The NASA STI Program

provides access to the NASA Aeronautics

and Space Database and its public

interface, the NASA Technical Report

Server, thus providing one of the largest

collection of aeronautical and space

science STI in the world. Results are

published in both non-NASA channels and

by NASA in the NASA STI Report Series,

which includes the following report types:

• TECHNICAL PUBLICATION. Reports

of completed research or a major

significant phase of research that present

the results of NASA programs and

include extensive data or theoretical

analysis. Includes compilations of

significant scientific and technical data

and information deemed to be of

continuing reference value. NASA

counterpart of peer-reviewed formal

professional papers, but having less

stringent limitations on manuscript

length and extent of graphic

presentations.

• TECHNICAL MEMORANDUM.

Scientific and technical findings that are

preliminary or of specialized interest,

e.g., quick release reports, working

papers, and bibliographies that contain

minimal annotation. Does not contain

extensive analysis.

• CONTRACTOR REPORT. Scientific

and technical findings by

NASA-sponsored contractors and

grantees.

• CONFERENCE PUBLICATION.

Collected papers from scientific and

technical conferences, symposia,

seminars, or other meetings sponsored

or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION.

English- language translations of foreign

scientific and technical material

pertinent to NASA’s mission.

Specialized services also include creating

custom thesauri, building customized

databases, and organizing and publishing

research results.

For more information about the NASA STI

Program, see the following:

• Access the NASA STI program home

page at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA STI

Help Desk at 443-757-5803

• Phone the NASA STI Help Desk at

443-757-5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace

Information

7115 Standard Drive

Hanover, MD 21076–1320
NASA/TM–2014–218314

Exploratory Study for

Continuous-Time Parameter

Estimation of Ankle Dynamics

Sunil L. Kukreja

NASA Armstrong Flight Research Center

Richard D. Boyle

NASA Ames Research Center

National Aeronautics and

Space Administration

NASA Armstrong Flight Research Center, Edwards, California 93524-0001

May 2014
The use of trademarks or names of manufacturers in this report is for accurate reporting and does not

constitute an offical endorsement, either expressed or implied, of such products or manufacturers by

the National Aeronautics and Space Administration.

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076-1320

443-757-5802
Abstract

Recently, a parallel pathway model to describe ankle dynamics was proposed.

This model provides a relationship between ankle angle and net ankle torque as

the sum of a linear and nonlinear contribution. A technique to identify parameters of this model in discrete-time has been developed. However, these parameters are a nonlinear combination of the continuous-time physiology, making

insight into the underlying physiology impossible. The stable and accurate estimation of continuous-time parameters is critical for accurate disease modeling,

clinical diagnosis, robotic control strategies, development of optimal exercise

protocols for longterm space exploration, sports medicine, etc.

This paper explores the development of a system identification technique

to estimate the continuous-time parameters of ankle dynamics. The effectiveness of this approach is assessed via simulation of a continuous-time model of

ankle dynamics with typical parameters found in clinical studies. The results

show that although this technique improves estimates, it does not provide robust estimates of continuous-time parameters of ankle dynamics. Due to this

we conclude that alternative modeling strategies and more advanced estimation

techniques be considered for future work.

1
Nomenclature

Acronyms

IRF- Impulse Response Function

NARMAX- Nonlinear AutoRegressive, Moving Average eXogenous

ELS - Extended Least-Squares

PEI - Prediction Error Identification

SLS - Standardized Least-Squares

PRBS - Pseudo-Random Binary Sequence

SNR - Signal-to-Noise Ratio

NF - Noise-Free

CT - Continuous-Time

FIR - Finite Impulse Response

IIR - Infinite Impulse Response

LTI - Linear Time Invariant

Symbols

Fl - nonlinear mapping

u - exogenous input

z - output

e - innovation, or uncontrolled input

n - sample index point

l - maximum polynomial order

nz - maximum output lag

nu - maximum input lag

ne - maximum error (innovation) lag

I- inertia

B - viscosity

K - elasticity

Δ - reflex time delay

ζ - damping parameter

ω - natural frequency

g - gain

c0 - zeroth-order coefficient of polynomial nonlinearity

c1 - first-order term coefficient of polynomial nonlinearity

c2 - second-order term coefficient of polynomial nonlinearity

s - Laplace variable

U(s) - ankle angle

V(s) - ankle velocity

V(s)e−Δs - delayed ankle velocity

X(s) - output of half-wave rectifier

YL(s) - noise-free linear path output

YNL(s) - noise-free nonlinear path output

Y(s) - noise-free net torque

2
E(s) - measurement noise

θ0 −θ11 - discrete-time coefficients of ankle dynamics model

τ - discrete-time delay

υ(n)- terms of discrete-time ankle model associated with nonlinear polynomial

coefficient c1

χ(n)- terms of discrete-time ankle model associated with nonlinear polynomial

coefficient c2

T - sampling rate

εˆ - prediction errors

R - set of real numbers

N - data length

Z - vector of measured outputs

Zˆ - vector of predicted outputs

θˆ - extended least-squares estimate of the system parameters

p - number of system parameters

Ψ - partitioned regressor matrix

Ψzu - regressor matrix that is a function of z and u only

Ψzuεˆ - regressor matrix that represents all the cross products involving εˆ

Ψεˆ - regressor matrix that is a polynomial function of the prediction errors

Ψ- centered and standardized variate of Ψ

Z- centered and standardized variate of Z

ΣΨ - diagonal matrix of standard deviations

ΣΨk - standard deviation of the kth column of Ψ

μΨ - matrix who’s kth column has all entries equal to the mean of column k of

Ψ

θ

- centered and standardized system parameters

minθ

1

2









Z−Ψθ









2

2

- standardized least-squares objective function

Hz - Hertz

ms - millisecond

dB - decibel

3
1 Introduction

In a zero or low-g environment it is critical to keep astronauts healthy and

functional. To maintain bone and muscle mass, astronauts need to exercise everyday [1]. Under terrestrial 1-g conditions gravity works against muscles and

bones, which requires the neuromuscular system to maintain enough muscle

and bone mass to support body weight. In an altered gravity-inertial environment the forces of gravity are significantly reduced. As a result, astronauts

lose muscle mass and bone density since it is not required to support body

weight [2].

Under current technological limitations (e.g.lack of artificial gravity through

inertial forces), rigorous physical activity is the only successful way to compensate for the lack of gravity [2]. Nevertheless, even with rigorous exercise, astronauts have typically lost 0.4-1% of their bone density per month in space [3].

Upon return to Earth, with a thorough physical therapy regime astronauts progressively recover muscle tissue and much of the bone mass lost during a semilong duration (i.e. 4-6 months) in low Earth orbit. However, throughout the

mission span it is important that astronauts are strong enough to perform strenuous activities in space, such as spacewalks and emergency procedures during

landing. Due to these possible vital mission activities, a regular exercise routine

in orbit prepares astronaut for such situations and accelerates a reconditioning

period to recover muscle and bone loss.

A parallel pathway model of ankle dynamics has been proposed that measures changes to intrinsic, reflex and muscle properties [4]. These characteristics change, for example, due to a lack of effective exercise. Therefore, astronauts health and exercise effectiveness can be monitored in orbit by evaluating

deviations of these parameters from optimal pre-orbit measures. In this work,

we propose to use this parallel pathway model as a foundation to assess lower

limb health during orbit to maintain an optimal exercise program.

Traditional approaches to nonlinear system identification of human ankle

dynamics have relied on quasi-linear methods, e.g. impulse response function

(IRF) method [5]. These methods provide convenient, robust means of characterizing the dynamics of nonlinear systems without requiring a priori assumptions regarding the system structure. However, nonparametric techniques may

require many parameters to describe even simple systems and can be difficult

to relate to the parameters of the underlying physiological system.

The NARMAX (Nonlinear AutoRegressive, Moving Average eXogenous)

model structure has been shown to be well suited to modeling the input-output

behavior of ankle dynamics [6]. The unknown model parameters can be estimated in discrete-time using the extended least-squares (ELS) algorithm [7–

9]. However, discrete-time parameters are a nonlinear combination of the

continuous-time physiology, making insight into the underlying mechanisms

difficult. The primary objective of biological system identification is to provide

parameters that give insight and relevance to the underlying biology, which are

represented in continuous-time.

The problem of continuous-time system identification from sampled input4
output data can be divided into two broad approaches: (i) indirect methods,

where a discrete-time model is estimated from sampled data; then an equivalent continuous-time model is calculated and (ii) direct methods, where a

continuous-time model is obtained directly without going through the intermediate step of first determining a discrete-time model; based on concepts

of approximate numerical integration to recreate time-derivatives needed in

continuous-time formulations [10].

The central problem for indirect methods is that of estimating stable and

robust continuous-time parameters from discrete-time estimates. For both linear and nonlinear systems, the inverse mapping problem often results in large

bias or sign flips. To address these issues several indirect techniques have

been developed to modify the estimators, keeping the computational complexity low [11–13]. In addition, a method for estimating the natural logarithm of

a square matrix to map discrete-time parameters to continuous-time has been

proposed [14]. Most other procedures are variants of these two approaches.

Direct methods avoid the inverse mapping problem by identifying a continuoustime model directly from sampled data. This raises several technical issues.

Unlike a difference equation model, a differential equation model contains

time-derivative terms that may be required but not available from measurements. Many techniques have been devised to deal with the need to reconstruct

these time-derivatives (see e.g. [10, 15–17]). However, these approaches are

challenged by a lack of robust implementation for numerical computation of

derivatives, selection of filter cutoffs to suppress high frequency noise, filter

order, etc. These problems are more complex when the system under study is

nonlinear, which is inherently the case for biological systems.

We propose to investigate the problem of continuous-time system identification by developing methodology using an indirect approach. One such approach relies on matrix preconditioning techniques, such as standardized leastsquares, to improve the spectral properties of the regressor matrix [18,19]. The

transformed matrix will have a smaller spectral condition number, and eigenvalues clustered around one. Often when the clustered spectrum is away from

zero it results in rapid and robust convergence, especially when the preconditioned matrix is close to normal. We deem this will provide more robust solutions and greater consistency for nonlinear biological systems by circumventing issues of implementing numerical derivatives and filter selection required

by direct techniques, which are more challenging in a nonlinear framework.

Here, we focus on one critical issue, namely, that of mapping the underlying

continuous-time system to discrete-time for estimation, then inverse mapping

back to continuous-time to provide physiological relevance and insight.

The organization of this paper is as follows. The NARMAX model structure is described in section 2. Section 3 discusses a continuous-time representation of a parallel pathway model describing ankle dynamics and its NARMAX

representation. In Section 4 we summarize the parameter estimation method

implemented in this study, namely, extended least-squares. Section 5 shows

how to standardized the identification data to formulate a standardized leastsquares approach, which we investigate for its utility as an indirect technique

5
to estimate continuous-time parameters of ankle dynamics. Section 6 illustrates

the results of our proposed approach on a simulated ankle model. Section 7 provides a discussion of our findings and suggestions for future direction. Section

8 summarizes the conclusions of our study.

2 NARMAX Model Description

The input-output relationship of many nonlinear dynamic systems can be written as the nonlinear difference equation

z(n) = Fl

[z(n−1),··· ,z(n−nz),u(n),··· ,u(n−nu), (1)

e(n−1),··· , e(n−ne)] +e(n)

where Fl is a nonlinear mapping, u is the exogenous input, z is the output,

and e is the innovation, or uncontrolled input. This structure describes both the

stochastic and deterministic components of a system. This nonlinear mapping

may include a variety of nonlinear terms, such as terms raised to an integer

power (e.g. u2(n−3)), products of present and past inputs (e.g. u(n)u(n−1)),

past outputs (e.g.z(n−1)z(n−2)), or cross-terms (e.g. u2(n−1)z(n−2)). This

system description encompasses most forms of nonlinear difference equations

that are linear-in-the-parameters.

3 Parallel Pathway Model of Ankle Dynamics

The Neuromuscular Control Laboratory (Department of Biomedical Engineering, McGill University, Montreal Canada), has developed a parallel pathway  ́

model (Fig. 1) to describe ankle dynamics [4]. This model provides a relationship between ankle angle (rad) and net ankle torque (Nm) as the sum of a

linear and nonlinear contribution. The upper, linear pathway models intrinsic

stiffness as a second-order system with parameters corresponding to inertia (I),

viscosity (B) and elasticity (K). The lower, nonlinear pathway models reflex

stiffness as a cascade of a derivative, a reflex time delay, a static nonlinearity

(i.e. half-wave rectifier) and a low-pass system representing muscle activation.

The parameters associated with the low-pass system are damping parameter

(ζ ), natural frequency (ω) and gain (g). This model is ideally suited to study

the effects of exercise (in space) because it quantifies both the intrinsic and reflex stiffness components. The parameters of these two paths will deviate from

normal a priori terrestrial measurements due to a suboptimal exercise regime.

A method to monitor changes of these parameters could be used to develop

individualized exercise protocols for astronauts.

Kukreja et al. [6] showed that a second-order static polynomial (c0+c1x(n)+

c2x2(n)) provided a good approximation to the half-wave rectifier which gave

6
✖✕

✗✔

Σ ✖✕

✗✔

Σ



✲ Δ ✲ V(s)e−Δs

❄

✲ ✲

❄

✻

✲

❄

✲ Is2 +Bs+K

+

Torque

Measured

Torque

Net +

+

+ Y(s)

Z(s)

E(s)

Reflex Delay

Static Nonlinearity

(Half-Wave Rectifier)

gω2

s2+2ζωs+ω2

Ankle Angle

U(s)

Ankle Velocity

V(s)

s

X(s)

Intrinsic Component

YL(s)

YNL(s)

Muscle Activation

Figure 1. Continuous-time morphological model describing ankle dynamics.

Upper (linear) path: intrinsic stiffness. Lower (nonlinear) path: reflex stiffness.

a NARMAX representation for this model of ankle dynamics as

z(n) = θ0 +θ1z(n−1) +θ2z(n−2) +θ3u(n) +θ4u(n−1) +θ5u(n−2) +θ6u(n−3(2) )

+ θ7u(n−4) +θ8[u(n−τ) +u(n−τ −1)−u(n−τ −2)−u(n−τ −3)]

+ θ9[u2

(n−τ) +3u2(n−τ −1) +3u2(n−τ −2) +u2(n−τ −3)

− 2u(n−τ)u(n−τ −1)−4u(n−τ −1)u(n−τ −2)−2u(n−τ −2)u(n−τ −3)]

+ θ10e(n−1) +θ11e(n−2) +e(n)

= θ0 +θ1y(n−1) +θ2y(n−2) +θ3u(n) +θ4u(n−1) +θ5u(n−2) +θ6u(n−3)

+ θ7u(n−4) +θ8υ(n) +θ9χ(n) +θ10e(n−1) +θ11e(n−2) +e(n)

where τ is the discrete-time delay,

υ(n) = u(n−τ) +u(n−τ −1)−u(n−τ −2)−u(n−τ −3) (3)

and

χ(n) = u2

(n−τ) +3u2(n−τ −1) +3u2(n−τ −2) +u2(n−τ −3) (4)

− 2u(n−τ)u(n−τ −1)−4u(n−τ −1)u(n−τ −2)

− 2u(n−τ −2)u(n−τ −3).

Table 1 gives the relationships of these discrete-time NARMAX parameters

(Eqn. 2) to the underlying continuous-time coefficients (the “extra” coefficient T denotes the sampling rate). An estimate of unknown system parameters

can be obtained using standard prediction error identification (PEI) techniques,

such as extended least-squares.

7
NARMAX Relationship to

Coefficient Continuous-time Coefficient

θ0

4c0gω2T2

4+ω2T2+4ζωT

θ1 − −8+2ω2T2

4+ω2T2+4ζωT

θ2 −−4ζωT+4+ω2T2

4+ω2T2+4ζωT

θ3

I

T2 + B

T +K

θ4 (−2I

T2 − B

T )−((− −8+2ω2T2

4+ω2T2+4ζωT )( IT2 + B

T +K))

θ5 ( I

T2 )−((− −8+2ω2T24+ω2T2+4ζωT )(−2IT2 − B

T ))−((−−4ζωT+4+ω2T2

4+ω2T2+4ζωT )( IT2 + B

T +K))

θ6 −(− −8+2ω2T2

4+ω2T2+4ζωT )( IT2 )−((−−4ζωT+4+ω2T24+ω2T2+4ζωT )(−2IT2 − B

T ))

θ7 −((−−4ζωT+4+ω2T2

4+ω2T2+4ζωT )( IT2 ))

θ8 gω2T2c1

(4+ω2T2+4ζωT)T

θ9 gω2T2c2

(4+ω2T2+4ζωT)T2

Table 1. Theoretical relationship of NARMAX model parameters to

continuous-time system coefficients.

4 Parameter Estimation

Many parameter estimation techniques are based on least-squares theory. However, ordinary least-squares algorithms for linear systems cannot be applied

here because they assume that the noise terms in the model are independent

and the regressor matrix is deterministic [20]. Both of these conditions are violated in Eqn. 2, when one considers that only a noisy measure of the output is

available.

Clearly, the regressors in Eqn. 2 contain correlated error terms and noisy

data. To obtain unbiased parameters other estimation techniques based on leastsquares must be used. One such method that is applicable to NARMAX model

estimation is extended least-squares. This technique provides unbiased parameter estimates and is iterative [19, 21, 22].

4.1 Extended Least-Squares

Extended least-squares is one method, appropriate for NARMAX models that

easily enables unbiased estimates to be computed. ELS addresses the bias problem by modeling the lagged errors to obtain unbiased parameter estimates.

ELS for linear systems has been widely studied and is also referred to as

Panuska’s method, the extended matrix method, or approximate maximum

likelihood [7–9].

In general, since the noise sequence is a realization of a stochastic process,

it is not possible to solve for the noise source e, and it will not be equal to the

8
prediction errors [23]. The prediction errors, εˆ ∈ RN×1, are defined as

εˆ = Z−Zˆ (5)

where Z ∈ RN×1 is the measured output and Zˆ = Ψθˆ ∈ RN×1 is the predicted

output. In ELS, the NARMAX formulation of Eqn. 1 is redefined into a prediction error model by replacing e with εˆ; making it a deterministic least-squares

problem.

The ELS formulation is defined as

θˆ = (ΨTΨ)

−1

ΨTZ; where Ψ = [ΨzuΨzuεˆΨεˆ]. (6)

θˆ ∈ Rp×1 is an ELS estimate of the system parameters, Ψ ∈ RN×p is a partitioned regressor matrix where Ψzu is a function of z and u only, Ψzuεˆ represents

all the cross products involving εˆ, and Ψεˆ is a polynomial function of the prediction errors only [23].

Often, ELS is considered a pseudolinear approach to parameter estimation

[7,9,24]. Strictly speaking, the introduction of prediction errors into the model

formulation no longer makes the model linear-in-the-parameters because the

prediction errors depend on the model output, which is a function of all model

parameters. The ELS technique solves a nonlinear optimization problem by

ignoring the nonlinear character of the model and employing a least-squares

approach. Essentially, ELS uses an approximate gradient of the model output

with respect to the model parameters as a regression vector.

5 Standardized Least-Squares

In many nonlinear systems there are large numerical differences of the regressors due to the basis function(s) used to estimate the static map. The

large differences lead to the regressor matrix being ill-conditioned and results

in unstable matrix inversion and poor parameter estimates [25]. To alleviate

ill-conditioning we propose using standardized least-squares (SLS). The SLS

technique is summarized as follows.

Given a matrix of independent variables Ψ and of dependent variables Z

compute the mean and standard deviation of each variable and replace Ψ and

Z with the centered and standardized variate as

Ψ= (Ψ− μΨ)Σ−1

Ψ and Z= (Z− μZ)Σ−1

Z (7)

where ΣΨ is a diagonal matrix of standard deviations with ΣΨk denoting the

standard deviation of the kth column of Ψ and μΨ is a matrix whose kth column

has all entries equal to the mean of column k of Ψ. Substituting Eqn. 7 into 6

yields a SLS formulation. The SLS objective function is compactly expressed

as

min

θ



1

2









Z−Ψθ









2

2

. (8)

This is an unbiased estimator and converges asymptotically to the true system

parameters.

9
6 Simulations

The efficacy of the SLS ELS technique was assessed to provide stable estimates

of continuous-time parameters describing ankle dynamics by simulating the

parallel pathway model in continuous-time using Simulink (Fig.1). The inputs

were bandlimited 30 Hz pseudo-random binary sequences (PRBS) which were

0.02 rad (peak-to-peak) and had a 35 ms switching rate. A PRBS input was

used to excite the system dynamics since it is standard practice in experimental

settings. The continuous-time coefficients used in this study correspond to

values found in experiments and are shown in Table 2. One thousand MonteCT Coefficient Value

I 0.015 Nm/rad/s2

B 0.8 Nm/rad/s

K 150 Nm/rad

ω 40.0

ζ 1.00

g 10.00 Nm/rad/s

Δ 0.045 s

T 0.005 s

Table 2. Continuous-time coefficient values. I: inertia, B: viscosity, K: elasticity, ω: natural frequency, ζ : damping parameter, g: reflex stiffness gain, Δ:

reflex delay and T: sampling interval.

Carlo simulations were generated in which each input-output realization was

unique, and had a unique Gaussian, white, zero-mean, noise sequence added to

the output. Excluding the noise-free (NF) case the signal-to-noise ratio (SNR)

of the noise sequence was decreased from 20 to 0 dB in increments of 5 dB.

For identification, the data length was N = 4,000 points.

6.1 Continuous-Time Parameter Estimation

The system parameters were estimated as outlined in Eqns. 6 – 8. Specifically,

for each input-output realization, we analyzed the ankle dynamics model as

follows:

1. Non-Standardized Extended Least-Squares: Ankle dynamics were analyzed using the ELS estimator.

2. Standardized Extended Least-Squares: Ankle dynamics were analyzed

using the ELS estimator and standardized data.

Fig. 2 shows a typical input-output realization used for this analysis. Employing

our approach, continuous-time parameters were computed from discrete-time

NARMAX estimates using the theoretical relationships given in Table 3.

10
Figure 2. Typical PRBS ankle position input and torque output with 20 dB

SNR used for Monte-Carlo analysis.

CT Coefficient DT Relationship

I = θ7 × T2

−θ2

B = θ6+θ1 I

T2

−θ2 = −( 2I

T2 + B

T ) + 2I

T2 = −B

T × −T

K = θ3 −( I

T2 + B

T )

ω = −4+4θ2+4θ1

−1+θ2−θ1 =



(ω2T2)× 1

T2

ζ = −2−2θ2

−1+θ2−θ1 = ζωTωT

g = θ8 ×(4+ω2T2 +4ζωT)×T = gc1ω2T2

c1ω2T2

Table 3. Discrete to continuous-time relationships for parameters I, B, K, ω, ζ

and g of the parallel pathway ankle model.

11
The results of this study are presented in Fig. 3. The panels show plots of

standard deviation about the mean for estimated continuous-time parameters of

the linear and nonlinear paths (I, B, K, g, ω and ζ ) using the non-standardized

and standardized ELS identification techniques. Column one shows the linear

path parameters. The continuous-time estimates using non-standardized ELS

and relationships given in Table 3 provided a mean value that was close to

the true parameter and had monotonically increasing variance from 20-10 dB

SNR. However, with decreasing SNR (5-0 dB) the parameter mean was biased

away from the true with non-monotonically increasing variance that included

estimates with incorrect sign. Using standardized data with the ELS algorithm,

mean and variance estimates for the linear path parameters improved for all

SNR levels but experienced a similar non-monotonically increasing variance

for 5-0 dB SNR and encompassed values with opposite sign. Column two

represents the nonlinear path estimates. Applying non-standardized and standardized ELS to the nonlinear path illustrates that both approaches gave similar results. For the damping (ζ ) and gain (g) parameters standardized data

provided estimate that were worse than non-standardized data. The natural frequency (ω) estimates were moderately improved using standardized data. Results show that, for the linear path, the parameters were close to the true mean

but variance estimate was non-monotonic and included estimates with incorrect sign for low SNR levels. However, for the nonlinear path, continuous-time

parameters computed using both approaches were significantly biased from the

true mean for all SNR levels but did not include parameter estimates with incorrect sign.

7 Discussion

In this section, we summarize the findings of this study and discuss possible

future approaches. The discussion highlights major features of our approach

and offers possible explanations for the unexpected errors. The suggestion for

future direction provided two broad alternatives to our solution, namely, an

alternative identification strategy and discretization approach.

7.1 Continuous-Time Parameter Estimation

Estimation of continuous-time parameters of ankle dynamics demonstrates that,

with noise-free input and output the both the ELS and SLS ELS approaches

provided estimates that agree with the true known values. However, with output additive noise the linear path parameter estimates had non-monotonically

increasing variance and encompassed values with opposite sign (see column 1,

Fig. 3). For the nonlinear path, the parameters were biased and the standard

deviation did not encompass the true value (see column 2, Fig. 3). Although

the nonlinear path parameters were biased, they did not experience sign flips as

the linear path parameters did.

The reason for the parameter bias may be due to the opposing nature of

12
Figure 3. Non-standardized and standardized ELS approaches to continuoustime parameter estimation of ankle dynamics. Column 1: Linear path parameters, I,B,K. Column 2: Nonlinear path parameters, ζ ,ω,g. Ordinate: STD

about mean. Abscissa: Output SNR = NF, 20, 15, 10, 5, 0 dB, where NF

denotes noise-free case. (Note that the abscissa is shown in decreasing SNR

which corresponds to increasing noise amplitude.)

13
the two paths. The linear path is described by a finite impulse response (FIR)

system, which has theoretical poles close to the Nyquist frequency; requiring a

high frequency signal to excite the dynamics sufficiently. The nonlinear path is

described by a combination of a nonlinear static function and linear dynamics.

The nonlinear function is described as a half-wave rectifier. This is considered a

hard nonlinearity because it generates multiple harmonics and its output cannot

be written analytically. To limit the effect of higher order harmonics and avoid

internal aliasing at the nonlinearity output the input signal was bandlimited to

be of low frequency. The opposing requirements of the two paths are central to

this estimation problem and difficult to satisfy simultaneously. It is not possible

to make more definitive comments about this approach without further detailed

investigation and analysis. Below we outline a few alternative approaches to

solve this identification problem.

7.2 Closed-Loop Simulation

A possible solution to continuous-time parameter estimation for ankle dynamics is to identify the model in closed-loop, removing the effects of derivatives

(FIR realization of the linear path). In feedback, intrinsic stiffness is represented as a compliance model (pure infinite impulse response (IIR) realization).

Although it may be possible to pose the identification problem as a feedback

model, the true system is represented from position input to torque output. The

recorded output torque has observation noise associated with it. When identifying the model in feedback the noisy output is considered the input. Noise

on the input violates assumptions and conditions for least-squares to yield unbiased parameter estimates. Implementing a feedback identification approach

may require advanced methods such as total-least squares [26, 27].

7.3 Discrete to Continuous-Time Parameter Mapping

Another possibility for the source of error for continuous-time parameters computed using our NARMAX approach may be related to the nonlinear relationships between discrete-time NARMAX parameters and continuous-time parameters of the physical system (see Table 3). A small deviation from the true

parameter value in discrete-time, due to noise or numerical error, may appear

as a significant error in the estimated continuous-time coefficients. As a result

it may be advantageous to study these model parameters only in discrete-time

as a combination of physiological effects.

7.4 Continuous to Discrete-Time Transformations

To derive a NARMAX model formulation the bilinear transform and Newton’s

backwards formula were used to convert the continuous-time linear dynamics to discrete-time. These transforms were implemented since both require

only a simple substitution to convert a continuous-time system to discrete-time.

Two other techniques that give better approximation for linear time invariant

14
(LTI) systems are linear extrapolation and linear interpolation methods [28,29].

However, these techniques are seldom used due to added complexity for little

gain.

Linear extrapolation gives an improper transfer function (i.e. more zeros

than poles) and linear interpolation gives a transfer function that is strictly

proper (i.e. equal zeros and poles) [28, 29]. These methods produce a discretetime transfer function that gives a better output response than the bilinear transform or Newton’s backwards formula. However, the pitfall of these methods

is that it is difficult to derive the coefficients of the discrete-time linear system.

This is the main reason that almost all engineering text books and literature

only discussed the bilinear transform and Newton’s backwards formula.

Using linear interpolation or linear extrapolation to derive a NARMAX

representation of ankle dynamics it may be possible to derive a better approximation to the derivative than the one used i.e. Newton’s backwards formula.

This may give a better approximation to the derivative and provide simulation

data that matches the continuous-time process better. However, one drawback

is that it may give a discrete-time approximation to the derivative that is higher

than first order thereby increasing the complexity of the NARMAX representation [28, 29]. Another limitation is that the continuous to discrete mapping

of the continuous-time parameters will be more complicated since it involves

exponential functions. This may result in more sensitivity for continuous-time

parameter estimation since small deviations from the true value in discrete-time

will result in large errors in continuous-time.

8 Conclusions

Clearly, much work remains to be done to resolve the possible source(s) of

error for continuous-time parameter estimation for the model of ankle dynamics. Our results show that the proposed standardized extended least-squares

method provided some improvement for estimating the continuous-time parameters of the linear path. However, the parameter variance was too high

because it included parameters with incorrect sign. For the nonlinear path, this

technique did not offer any improvement over the non-standardized extended

least-squares approach. A possible explanation for these results is due to opposing constrains needed to properly excite the dynamics of both paths simultaneously. The linear path requires a high frequency signal while the nonlinear

path needs low frequency. It is unclear how to satisfy both concurrently. It

is not possible to make more definitive remarks about the root of this error

and remains an open research question. We have offered several suggestions

for alternative ways forward to estimate the continuous-time parameters of the

parallel pathway model of ankle dynamics. Of these, we deem that identifying

this model in closed loop and implementing more advanced estimation techniques such as total least-squares is likely the best way forward.

15
References

1. V. Convertino, “Exercise as a countermeasure for physiological adaptation to prolonged spaceflight,” Medicine & Science in Sports & Exercise,

vol. 28, no. 8, pp. 999 –1014., 1972.

2. K. Baldwin, T. White, S. Arnaud, V. Edgerton, W. Kraemer, R. Kram,

D. Raab-Cullen, and C. Snow, “Musculoskeletal adaptations to weightlessness and development of effective countermeasures,” Medicine & Science

in Sports & Exercise, vol. 28, no. 10, pp. 1247–1253, 1996.

3. J. Iwamoto, T. Takeda, and Y. Sato, “Interventions to prevent bone loss in

astronauts during space flight,” Keio Journal of Medicine, vol. 54, no. 2,

pp. 55–59, 2005.

4. R. Kearney, R. Stein, and L. Parameswaran, “Identification of intrinsic

and reflex contributions to human ankle stiffness dynamics,” IEEE Trans.

Biomed. Eng., vol. 44, no. 6, pp. 493–504, 1997.

5. R. Kearney and I. Hunter, “System identification of human joint dynamics,” Critical Reviews in Biomedical Engineering, vol. 18, no. 1, pp. 55–87,

1990.

6. S. Kukreja, H. Galiana, and R. Kearney, “NARMAX representation and

identification of ankle dynamics,” IEEE-TBME, vol. 50, no. 1, pp. 70–81,

2003.

7. V. Panuska, “A stochastic approximation method for identification of linear systems using adaptive filtering,” in Proceedings 9th Joint Automatic

Control Conference. Ann Arbor, MI, USA: IEEE, June 1968, pp. 1014–

1021.

8. ——, “An adaptive recursive least squares identification algorithm,” in

Proceedings 8th IEEE Symp. on Adaptive Processes, University Park, PA,

USA, Nov. 1969, p. paper 6e.

9. P. Young, “The use of linear regression and relaxed procedures for the

identification of dynamic processes,” in Proceedings 7th IEEE Symp. on

Adaptive Processes, Los Angeles, CA, USA, Dec. 1968, pp. 501–505.

10. E. D. Sontag, in Identification of continuous-time systems. Methodology

and computer implementation, N. Sinha and G. Rao, Eds. Kluwer Academic Publishers, 1991.

11. E. K. Larsson and T. Soderstr  ̈ om, “Continuous-time ar parameter estima-  ̈

tion by using properties of sampled systems,” in Proceedings 15th IFAC

World Congress, vol. 15, Barcelona, Spain, July 2002.

12. H. Fan, T. Soderstr  ̈ om, M. Mossberg, B. Carlsson, and Y.Zou, “Estimation  ̈

of continuous-time ar process parameters from discrete-time data,” IEEE

Transactions on Signal Processing, vol. 47, no. 5, pp. 1232–1244, 1999.

16
13. T. Soderstr  ̈ om, H. Fan, B. Carlsson, and S. Bigi, “Least squares parameter  ̈

estimation of continuous-time arx models from discrete-time data,” IEEE

Transactions on Automatic Control, vol. 42, no. 5, pp. 659–673, 1997.

14. G.J.Lastman and N. Sinha, “Continuous-time model identification from

sampled data. implementation issues and performance evaluation,” Electronics Letters, vol. 27, no. 16, pp. 1468–1469, 1991.

15. H. Garnier, M. Mensler, and A. Richard, “Infinite series for logarithm of

a matrix, applied to identification of continuous-time multivariable systems from discrete-time models,” International Journal of Control, vol. 76,

no. 13, pp. 1337–1357, 2003.

16. P. Young, H. Garnier, and A. Jarvis, “The identification of continuous-time

linear and nonlinear models: a tutorial with environmental applications,”

in 13th IFAC Symposium on System Identification (SYSID2003)s, vol. 13,

Rotterdam, The Netherlands, August 2003, pp. 618–628.

17. Unbehauen and G. Rao, “A review of identification in continuous-time

systems,” Annual reviews in Control, vol. 22, pp. 145–171, 1998.

18. L. N. Trefethen and D. Bau, Numerical Linear Algebra. Philadelphia,

PA: SIAM, 1997.

19. L. Ljung, System Identification: Theory for the User, 2nd ed. Upper

Saddle River, New Jersey: Prentice Hall, Inc., 1999.

20. T. C. Hsia, System Identification: Least-Squares Methods. Lexington,

Massachusetts: Lexington Books, D.C. Heath and Company, 1977.

21. S. Chen and S. Billings, “Representations of non-linear systems: the NARMAX model,” Int. J. Control, vol. 49, no. 3, pp. 1013–1032, 1989.

22. T. Soderstr  ̈ om and P. Stoica,  ̈ System Identification. Prentice-Hall International, 1989.

23. S. Billings and W. Voon, “Least squares parameter estimation algorithms

for non-linear systems,” Int.J. System Sci., vol. 15, no. 6, pp. 601–615,

1984.

24. G. Goodwin and R. Payne, Dynamic System Identification: Experiment

Design and Data Analysis, ser. Mathematics in Science and Engineering.

New York: Academic Press, 1977, vol. 136.

25. P. Hansen, V. Pereyra, and G. Scherer, Least Squares Data Fitting with

Applications, 1st ed. Baltimore, MD: John Hopkins University Press,

2012.

26. A. Bj  ̊ orck,  ̈ Numerical Methods for Least Squares Problems, 1st ed.

Philadelphia: Siam, 1996.

17
27. S. V. Huffel and J. Vandewalle, The Total Least Squares Problem: Computational Aspects and Analysis, 1st ed. Philadelphia: Siam, 1991.

28. S. Haykin, “A unified treatment of recursive digital filtering,” IEEE Trans.

Auto. Control, vol. 17, pp. 113–116, 1972.

29. N. K. Sinha, “Estimation of transfer function of continuous system from

sampled data,” Proceeding of the IEE, vol. 119, no. 5, pp. 612–614, 1972.

18

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704–0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection

of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports

(0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be

subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

1. REPORT DATE (DD-MM-YYYY)

01-05-2014

2. REPORT TYPE

Technical Memorandum

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Exploratory Study for Continuous-Time Parameter Estimation of Ankle Dynamics

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

6. AUTHOR(S)

Sunil L. Kukreja and Richard D. Boyle

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA NASA Armstrong Flight Research Center, Edwards, California 93524-0001

8. PERFORMING ORGANIZATION

REPORT NUMBER

L–

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

10. SPONSOR/MONITOR’S ACRONYM(S)

NASA

11. SPONSOR/MONITOR’S REPORT

NUMBER(S)

NASA/TM–2014–218314

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category

Availability: NASA CASI (443) 757-5802

13. SUPPLEMENTARY NOTES

An electronic version can be found at http://ntrs.nasa.gov.

14. ABSTRACT

Recently, a parallel pathway model to describe ankle dynamics was proposed. This model provides a relationship between ankle angle and

net ankle torque as the sum of a linear and nonlinear contribution. A technique to identify parameters of this model in discrete-time has been

developed. However, these parameters are a nonlinear combination of the continuous-time physiology, making insight into the underlying

physiology impossible. The stable and accurate estimation of continuous-time parameters is critical for accurate disease modeling, clinical

diagnosis, robotic control strategies, development of optimal exercise protocols for longterm space exploration, sports medicine, etc.

This paper explores the development of a system identification technique to estimate the continuous-time parameters of ankle dynamics. The

effectiveness of this approach is assessed via simulation of a continuous-time model of ankle dynamics with typical parameters found in

clinical studies. The results show that although this technique improves estimates, it does not provide robust estimates of continuous-time

parameters of ankle dynamics. Due to this we conclude that alternative modeling strategies and more advanced estimation techniques be

considered for future work.

15. SUBJECT TERMS

System Identification, Indirect System, Least-Squares, Standardized Least-Squares, Switched Systems, Parametric Models, NARMAX,

Exercise Protocols

16. SECURITY CLASSIFICATION OF:

a. REPORT

U

b. ABSTRACT

U

c. THIS PAGE

U

17. LIMITATION OF

ABSTRACT

UU

18. NUMBER

OF

PAGES

16

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802