Convex Computation of the Basin of Stability to
Measure the Likelihood of Falling: A Case Study
on the Sit-to-Stand Task
Victor Shia∗ , Talia Moore† , Ruzena Bajcsy∗ , Ram Vasudevan‡

arXiv:1604.00602v1 [] 3 Apr 2016

∗ Electrical

Engineering and Computer Sciences
University of California, Berkeley
{vshia, bajcsy}@berkeley.edu
† Organismic and Evolutionary Biology
Harvard University
taliaym@gmail.com
† Mechanical Engineering
University of Michigan
ramv@umich.edu

Abstract—Locomotion in the real world involves unexpected
perturbations, and therefore requires strategies to maintain
stability to successfully execute desired behaviours. Ensuring the
safety of locomoting systems therefore necessitates a quantitative
metric for stability. Due to the difficulty of determining the set
of perturbations that induce failure, researchers have used a
variety of features as a proxy to describe stability. This paper
utilises recent advances in dynamical systems theory to develop
a personalised, automated framework to compute the set of
perturbations from which a system can avoid failure, which is
known as the basin of stability. The approach tracks human
motion to synthesise a control input that is analysed to measure
the basin of stability. The utility of this analysis is verified on a
Sit-to-Stand task performed by 15 individuals. The experiment
illustrates that the computed basin of stability for each individual
can successfully differentiate between less and more stable Sitto-Stand strategies.
Index Terms—Stability analysis, Locomotion Biomechanics,
Optimization and Optimal Control, Sit-to-Stand

I. I NTRODUCTION
Falls are a leading cause of accidental injury and death
throughout much of the world. Due to the aging population
and the outsized impact falling has on the elderly, the cost
associated with falls is expected to rise dramatically in the next
twenty years [1]. Directed therapeutic care can significantly
reduce the risk of falling [2, 3]; however, the resources
available for such treatment are limited. An automated test
identifying individuals at risk for falling can make targeted
deployment of therapeutic care feasible. Unfortunately the
construction of such a test has been challenging.
This paper develops a personalised automated diagnostic
test that uses kinematic observations to measure an individual’s
likelihood of falling. The approach, which is grounded in
dynamical systems theory, computes the Basin of Stability
(BOS) of a locomotor pattern, or the set of perturbations that
do not lead to a fall under an individual’s chosen locomotor
strategy (illustrated in Figure 1). Informally, an individual that

is able to tolerate a larger set of perturbations has a larger BOS
and is less likely to fall.
In fact, measuring the BOS is a direct way to characterise
the likelihood of falling, since it identifies the specific deficiencies that lead to failure [4]. Unfortunately the computation
of this individual- and behaviour-specific BOS is challenging,
since it requires measuring the effect of arbitrary perturbations
to a nonlinear system. An empirical experiment would require
exhaustive perturbation of an individual throughout a locomotor pattern, which is practically infeasible and dangerous.
To address these issues, the presented approach computes
the BOS in a tractable manner using convex optimization.
Though the method is applicable to arbitrary locomotor patterns, this paper illustrates the utility of this technique by
analysing Sit-to-Stand (STS) manoeuvres, STS manoeuvres
are less complex than other locomotor patterns (e.g. walking,
running, climbing, lifting), simplifying the validation of the
method. Although comparatively straightforward, the ability
to stand is a prerequisite for bathing, cooking, dressing,
maintaining hygiene, and walking. As a result, difficulty in
performing STS manoeuvres is considered a primary risk
factor for falls amongst the elderly [5].
A. Existing Stability Metrics
Due to the importance of STS manoeuvres in maintaining
quality of life, and the impossibility of testing all possible
perturbations, a variety of methods to characterise an individual’s likelihood of falling while performing STS have been
proposed. These methods are summarised in Table I.
These methods generally summarise STS motions using a
single feature and perform versions of regression analysis to
estimate a patient’s stability. In doing so, they forfeit the ability
1 Note

that at time T , the BOS of the controller is small because we are
concerned about a finite-time BOS.

Methods
BERG Balance Test
Stops Walking When Talking
Timed Up & Go
Model based methods

Summary
A battery of functional tests with a single number that determines the likelihood of falling
Relates the amount of attention a person requires to perform an action with a likelihood of falling
Correlates a likelihood of falling with the the amount of time it takes to stand up
Uses a single inverted pendulum to determine the set of feasible initial positions or positions and
velocities that can lead to standing up

References
[6]
[7]
[8, 9]
[10, 11, 12,
13, 14, 15]

TABLE I: Various STS Stability Methods

Algorithm 1 Computing the Reachable Set for STS
(a) Pictorial Model

1:
2:
3:
4:
5:

(b) Basin of Stability over Time
Fig. 1: The BOS of an inverted pendulum swing-up controller with perturbations indicated by arrows. (a) illustrates a nominal (black) and perturbed
(magenta and cyan) swing-up trajectory for the inverted pendulum. (b)
illustrates the BOS (gray) of the nominal controller, the nominal trajectory
(black), and the perturbed trajectory (magenta and cyan) in the configuration
space of the inverted pendulum1 . Despite perturbation, the magenta and cyan
trajectories arrive at the upright configuration because they remain within the
basin of stability of the nominal trajectory.

to characterise the specific deficiencies limiting an individual’s
STS ability. More troublingly, according to several studies, the
ability of these clinical tests to distinguish between stable and
unstable patients is unclear [7, 16, 17, 18, 19].
B. Dynamical Systems Perspective
Computing the BOS is a fundamental objective of the
dynamical systems community since it can be used to verify
the satisfactory operation of a system despite perturbations.
In particular, engineers have long sought to understand the
behavior of a dynamical system after arbitrary perturbation.
Although it is possible to directly simulate arbitrary configurations (analogous to exhaustive perturbation), this technique
provides only limited insight. To address this issue, the dynamical systems community has studied numerical methods
to compute the BOS.
These methods for nonlinear systems include Lyapunovbased techniques [20, 21] and Hamilton-Jacobi based methods [22]. Lyapunov-based methods [23] search for functions
whose sub-level sets satisfy certain criteria. The construction
of such a Lyapunov function is possible for polynomial
dynamical systems using semidefinite programming [24], but
requires solving a challenging bilinear optimization problem, limiting potential applicability. Hamilton-Jacobi based
methods discretise the domain and run variants of dynamic
programming on a discretised nonlinear partial differential

Given: observations xobs of motion
Choose a model for the STS motion (Section II-B).
Run optimal control to find uobs (Section II-C).
Construct controller to track xobs (Section II-D).
Compute the backwards reachable set (Section II-E).

equation to determine the set of states that belong to a BOS.
Though this method is able to tractably compute the BOS
for dynamical systems with special structure [25], it is only
able to accurately compute the BOS for general systems with
less than 4 states. Recently, the authors developed a method
to analytically compute the BOS for polynomial dynamical
systems based on occupation measures. This method, which
relies on convex optimization and is described in further detail
below, tractably outer approximates the BOS of a system
without relying upon exhaustive perturbative experiments or
simulation. Furthermore, this method successfully synthesises
safe robotic motion for systems with up to 8 states [26, 27,
28, 29].
C. Contribution
First, Section II describes a personalised computational
framework to model, identify, and analyse the unique stability
of an individual’s motion. Second, Sections III and IV describe
a motion capture dataset of humans performing various STS
strategies. The proposed methods then evaluate each individual’s kinematic stability. Section V summarises the impact of
the proposed method and describes potential extensions.
II. M ETHODOLOGY
This section presents the framework to compute the BOS
of an individual’s locomotor pattern given kinematic data.
The approach is summarised informally in Algorithm 1.
The steps are described abstractly in this section to ensure
straightforward generalization to arbitrary locomotor patterns.
In Sections III and IV, a concrete instantiation of each step is
described in the case of STS motion.
A. Preliminaries
The notation used throughout the remainder of this paper
is presented in this section. Let Rn be a n-dimensional set
of real numbers. Let X ⊂ Rn be a compact set. Let [0, T ]
denote a time interval of interest. Let C 1 (X, R) be the space
of continuously differentiable functions from X to R. Let
L2 (X, R) be the space of square integrable functions under

the Lebesgue measure from X to R. Let Rn [x] be the set of
polynomials in x with maximum total degree n.
B. Model and Observations
Next, suppose that the dynamical model describing the
motion of an individual is:
ẋ(t) = fφ (t, x) + gφ (t, x)u(t, x)
x ∈ [x, x] ⊂ Rn
u ∈ [u, u] ⊂ R

(1)

m

where X = [x, x] ⊂ Rn represents the state space of the
model, f ∈ C 1 ([0, T ] × Rn , Rn ) and g ∈ C 1 ([0, T ] × Rn , Rm )
describe how the input u ∈ L2 ([0, T ] × Rn , Rm ) affect the
dynamics, φ represents the individual specific parameters of
the model (e.g. mass, limb length, moment of inertia, etc.),
and u, u ∈ Rm represent input bounds. As each individual is
different, φ, x, x, u and u are distinct for each individual and
must be identified as described in further detail in Section III.
This paper assumes that direct observations of the state trajectory of a nominal locomotor pattern, xobs ∈ C 1 ([0, T ], Rn ),
are available. This can be constructed after interpolation from
a variety of data sources as described in Section III.
C. Identifying an Input from Observations
After selecting a model, the input, uobs : [0, T ] → Rm that
generates the given observations must be constructed. There
are two methods for determining the input for the observed
motion: inverse dynamics and optimal control. Inverse dynamics uses the observed variables xobs : [0, T ] → X to estimate
ẋobs : [0, T ] → Rn . uobs can then be computed for all t in
Equation (1) using a known (xobs , ẋobs ) [30]. As the inverse
kinematic solution is sensitive to noise in xobs , optimal control
is used in this paper to compute uobs . Optimal control instead
calculates uobs via the optimization problem:
Z
inf
2

uobs ∈L ([0,T ],Rm )

T

kx(t) − xobs (t)k22 dt

(2)

0

s.t. ẋ(t) = fφ (t, x) + gφ (t, x)uobs (t)

∀t ∈ [0, T ]

x(t) ∈ X

∀t ∈ [0, T ]

uobs (t) ∈ [u, u]

∀t ∈ [0, T ]

The solution to this problem is a feedforward open loop
control input uobs that minimises the L2 error between the
state trajectory and the observed trajectory. After treating the
nominal input uobs as a polynomial function, collocation [31] is
used to transform this optimal control problem into a nonlinear
optimization program, which can be efficiently solved by a
variety of nonlinear programming solvers.
D. Feedback Controller Design
Neuroscientists, psychologists, motor control researchers,
and biomechanists have observed that the nominal trajectories
humans follow during locomotor patterns are robust to small
perturbations [32, 33, 34, 35, 36]. This robustness is conferred
by feedback about the nominal control input or goal. To date,

experimental research has been unable to identify an overall
strategy that endows such robustness.
For example, research has shown that subjects minimise
the square of jerk during reaching tasks [32]. Alternatively,
others have shown that for endpoint reaching tasks, subjects
utilise a time-varying Proportional Derivative (PD) control to
reach a specified endpoint [37]. For the lower body, other
researchers tracked the evolution of step width and found that
subjects tended to correct deviations with just a proportional
controller [38, 39].
To imbue the feedforward nominal control input that is
identified by the optimal control algorithm in Section II-C
with this feedback robustness, the following assumptions are
made:
Assumption 1. For each distinct locomotion action, humans
utilise a feedforward control law with corresponding feedback.
To perform a different action, the subject switches control laws.
Assumption 1 states that for a particular action, such as standing slowly, the subject follows a combination of feedforward
and feedback control laws. If a specific control law is not able
to take a subject to standing after perturbation, a subject must
switch control laws to stand safely.
Assumption 2. While performing a specific action, the subject
utilises a PD feedback around a nominal trajectory, xobs to
correct deviations in the trajectory.
According to Assumption 2 the feedback control law is:
u(t, x) = uobs (t) + ucc (t, x)

(3)

= uobs (t) + K(x(t) − xobs (t))
where ucc represents the general form of the feedback
controller and K represents the PD controller gain acting
on the states and observation. Note, the method presented
to estimate the BOS (described in Section II-E) can handle
more general nonlinear feedback control inputs. However, as
described earlier, the existing literature suggests that humans
apply only linear feedback [38, 39].
If the gain K on the feedback controller is selected too
rigidly, then the control law will oscillate around the desired
trajectory rather than converging to the final state of the desired
trajectory. This can be avoided with sufficiently small gains,
as illustrated by the system in Figure 1. To determine this
satisfactory feedback gain K, we apply a Linear Quadratic
Regulator (LQR) algorithm to determine the optimal state
feedback law u that minimises a quadratic cost:
Z T
min
(x(t) − xobs (t))T Q(x(t) − xobs (t))
ucc ∈L2 ([0,T ],Rm ) 0

+ ucc (t)T Rucc (t) dt (4)
s.t.

ẋ(t) = Ax(t) + B(uobs (t) + ucc (t, x)) ∀t ∈ [0, T ]

By selecting Q = I2 and R = 0.005I where I is the identity
matrix of appropriate dimension, the resulting controller is

designed to minimise the Q-weighted L2 error of x from
xobs . For linear systems, the LQR problem has a closed form
solution provided by the Algebraic Ricatti Equation described
by a linear state feedback law [40]. For the purposes of this
paper, small-angle approximations are used to linearise fφ and
gφ to obtain A and B.
Assumption 3. The torque limits are constant throughout the
motion.
As humans do not have the ability to apply arbitrary torque to
any joint, individual-specific torque limits [u, u] are set to the
minimum and maximum of uobs generated from the optimal
control.
E. Computing the Basin of Stability
Given a model, input bounds, and feedback control input
that tracks a nominal observation, the BOS can be formally
defined as follows: the BOS is the set of states as a function
of time that can be driven by the feedback control input to
a target configuration, XT ⊂ X, by time T . In the case of
STS, the target set XT corresponds to the set of states where
the subject is standing. For brevity, a modified optimization
algorithm inspired by [26] to compute this BOS is presented:

The intuition of the proof is as follows: if v(T, ·) ≥ α on XT
(5c), since v must decrease as the system evolves (5a), for a
point (t, x) ∈ [0, T ]×X to reach XT , v(t, x(t)) ≥ α must hold
for all time. As a result, the α super-level set of v at each time
t in [0, T ] can be used as a test to determine whether a point
does not belong to the BOS of the motion under consideration.
In Figure 1, for example, the light gray region denotes the
v(t, x) ≥ α level set with the dark gray region denoting
different time slices of the v(t, x) ≥ α level set. Outside
of the points that belong to BOS, the optimization problem
tries to minimise v by bringing it as close to 0 as possible.
Several recent papers formally describe the convergence of
this approach, which we do not include here for the sake of
brevity [26, 27, 28].
To solve (D) numerically, the dynamics are assumed to be
polynomial and the state space and target set are assumed to be
semi-algebraic sets. Since by the Stone-Weierstrass Theorem
polynomial functions are able to approximate the behavior
of other continuous functions on a compact domain [41],
this assumption is made without too much loss in generality.
The positivity constraints are converted to sum-of-squares
constraint [24]. The result is a semidefinite optimization program that tractably constructs an outer approximation to the
BOS [28].

Z
inf

v∈C 1 ([0,T ]×Rn ,Rn )
[0,T ]×X

s.t.

v(t, x)dtdx

(D)

III. E XPERIMENTS

∂v(t, x)
(fφ (t, x) + gφ (t, x)u(t, x))
(5a)
∂x
∂v(t, x)
+
≤ 0 ∀(t, x) ∈ [0, T ] × X
∂t
v(t, x) ≥ 0
∀(t, x) ∈ [0, T ] × X (5b)
v(T, x) ≥ α

∀x ∈ XT

(5c)

where α > 0 is a parameter that can be selected by the user.
To understand the relationship between the solution to this
optimization problem (v) and the BOS, notice that v(t, x) ≥ α
for points that belong on the BOS:
Lemma 1. If v is a solution to (D), then v(t, ·) ≥ α on the
BOS.
Proof. Suppose x : [0, T ] → X is a trajectory of the model
that reaches XT . Notice that x(t) is in the BOS for all t ∈
[0, T ]. Select an arbitrary ts ∈ [0, T ], then:
α ≤ v(T, x(T ))

(6)
Z

T



∂v(t, x(t))
(fφ (t, x(t)))
(7)
∂x
ts

∂v(t, x(t))
∂v
+
(gφ (t, x(t))u(t, x(t))) +
(t, x(t)) dt
∂x
∂t
(8)

= v(ts , x(ts )) +

≤ v(ts , x(ts ))

(9)

∂v(t,x(t))
since
(fφ (t, x(t)) + gφ (t, x(t))u(t, x(t))) +
∂x
∂v
∂t (t, x(t)) ≤ 0 on [0, T ] × X and v(T, ·) ≥ α on XT . The
desired result follows.

This section describes a formal implementation of the
method presented in Section II and an experiment constructed
to evaluate its validity. One method to verify the correctness of
a computed BOS is via direct perturbative experiments; however, these experiments can be prohibitive and are dangerous.
Instead we utilise observations from motor control research
to validate the computed stability estimates of distinct STS
maneuvers performed by each subject.
A. Intuition from Motor Control
Due to the time delay of the nervous system, motor
control researchers have hypothesised that the response of
perturbations to fast motions is largely governed by openloop reflex responses [42, 43], whereas slower motions allow
a closed-loop correcting response to perturbations. Based
on this experimentally validated tradeoff between speed and
feedback [44, 45, 46], we expect slower movements to have a
larger basin of stability.
The open- and closed-loop control laws are exemplified by
two distinct STS strategies: momentum-transfer and quasistatic [47, 48], shown in Figure 2. The momentum-transfer
strategy (indicated throughout in orange) consists of swinging
one’s trunk forward rapidly, using the forward momentum of
the upper body to stand up. This strategy requires significant
postural control, due to a dynamically unstable transition phase
[47]. In contrast, the quasi-static strategy (indicated throughout
in green) consists of leaning forward while sitting to position
the centre of mass (COM) above the feet, then using as
little momentum as possible to slowly stand. The motion
is statically stable at any given moment, but requires more

energy to perform than the momentum-transfer strategy [49].
Natural STS movements likely form a continuum between the
open-loop momentum transfer and the closed-loop quasi-static
strategies.
To validate Algorithm 1 experimentally, subjects performed
STS using their preferred strategy at two speeds and the momentum transfer and quasi-static strategies. Computed results
are considered accurate if (1) the slower preferred strategy
has a larger BOS than the faster preferred strategy and (2)
the quasi-static strategy has a larger BOS than the dynamic
strategy for the same individual.

(a) Seated subject prior to STS motion with (b) Circles indicate markers
feet on the force plate.
placed on front, Xs indicate
markers placed on back.

5

Knee
Hip

θ̇ (rad/s)

2.5

Fig. 3: Experimental Setup with LED Marker Placement

0

-2.5

-5
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Group

Time (s)

(a) Momentum Transfer STS diagram (top) and velocity profile (bottom).

Young

5

Knee
Hip

θ̇ (rad/s)

2.5

0

-2.5

Older

-5
0

0.5

1

1.5

2

2.5

3

Time (s)

(b) Quasi-static STS diagram (top) and velocity profile (bottom).

Subjects began in a seated position with their trunk and
tibiae oriented vertically, and arms crossed. The chair height
was adjusted such that the subject’s femurs were parallel to the
ground. Subjects wore a customised motion capture suit with
43 PhaseSpace markers (shown in Figure 3b). STS movements
were recorded using an AMTI OPT464508 force plate [50]
under the subject’s feet2 and a PhaseSpace Impulse X2 motion
capture system with 8 infrared cameras [51] (Figure 3a). Force
data were collected at 2400Hz, motion capture data were
collected at 480Hz, and the subject’s skeleton was extracted
using PhaseSpace’s Recap2 software [52]. Both the motion
capture and force plate data were smoothed using a 4th-order
Butterworth filter with a cut-off frequency of 2Hz.

measurements were used to determine the start and end time of the
STS motion

Age
23
25
37
26
22
28
29
21
25
25
84
69
77
76
74

Height (cm)
153.7
165.1
184.2
180.3
165.1
175.3
175.3
167.6
172.7
160.0
162.6
185.4
170.2
164.7
175.3

Weight (kg)
70.3
68.6
74.0
66.2
58.7
55.1
79.8
64.9
69.1
54.2
65.3
92.5
66.2
64.6
69.1

We collected data from 2 cohorts: 10 young and healthy
subjects and 5 older and healthy subjects3 . Data for individual
subjects is shown in Table II. Initially, subjects were asked to
stand without instruction to record the natural STS strategies
at slow and fast speeds (‘Untrained’ dataset). Subsequently,
subjects were shown videos demonstrating the momentum
transfer and quasi-static STS strategies to avoid individual
interpretation of the motion. The subjects were then asked to
perform the momentum transfer and quasi-static STS strategies
in a randomised order (‘Trained’ dataset’)4 .
C. Standing Models
The inverted pendulum model (IPM), shown in Figure 4a,
and the double inverted pendulum (DPM), shown in Figure 4b
are the controlled dynamic models for STS investigated in this
paper. Although the DPM is a more accurate representation of

3 This
2 Force

Gender
F
F
M
M
F
M
M
M
M
F
F
M
M
F
M

TABLE II: Data for each subject

Fig. 2: An illustrate of the two STS Strategies used to perform validation.
The difference in angular velocity of the ankle joint is negligible between the
two motions and is not shown.

B. Data Collection

ID
1
2
4
5
6
7
8
9
10
14
3
11
12
13
15

study was approved by the UC Berkeley Center for Protection of Human
Subjects, Protocol #2015-07-7767 and informed consent was obtained from
all subjects.
4 Videos can be found at: https://www.w3id.org/people/vshia/jrsi

θ2

r2

τ2

m1

l

θ

θ1

τ
l1

r1

τ1

(a) Inverted Pendulum Model (IPM) (b) Double Pendulum Model (DPM)
Fig. 4: The pair of STS controlled dynamic models considered in this paper.

human morphology, IPM has more widespread use due to the
complexity of DPM.
The IPM consists of an inverted pendulum attached to a
fixed foot on the ground with the point mass m at length l
away from the joint. Let θ represent the angle (with respect to
the vertical) and θ̇ represent angular velocity of the pendulum.
Together both variables are the state space of the IPM, while
τ represents the actuation at the ankle. The dynamics of the
IPM can be found in [10]. A 5th order Taylor expansion of
the dynamics is used while solving for the BOS using (D).
The motion capture of each individual was fit to the IPM by
setting m as the subject’s mass, l as the average distance from
the subject’s ankle to the subject’s COM, and θ as the angle
from the ankle to the subject’s COM.
The DPM is a double inverted pendulum attached to a fixed
foot on the ground. Let θ1 and θ̇1 represent the angle and
angular velocity of the lower link to the vertical and θ2 and θ̇2
represent the angle and angular velocity of the lower to upper
link shown in Figure 4b. Together these variables represent the
state space of the DPM while τ1 and τ2 represent the ankle
and hip actuation, respectively. The dynamics of the DPM can
be found in [53].
The motion capture of each individual was fit to the DPM
by setting m1 as the mass of the subject’s lower body (calf
and thigh), m2 as the mass of the subject’s upper body, l1
as the average length of the subject’s ankle to hip, r1 as the
average length from the ankle to the COM of the lower body,
and r2 as the average length from the hip to the COM of the
upper body. θ1 represents the angle from the subject’s ankle
to hip and θ2 represents the angle of the subject’s hip to upper
body.
For both models, masses and the COM positions of each
individual limb were computed using tabulated values found
in [54]. Individualised torque bounds are set to the maximum
and minimum torques obtained via the optimal control. The
domain bounds are the minimum and maximum observed
values from the data.

A. Results using an existing stability metric
The existing model-based approach for determining the
BOS of STS is called Region of Stability based on Velocity (ROSv). This method plots the normalised position and
velocity5 of the subject’s COM at the instant they rise from
the chair (as in Figure 5) [14, 15]. Points left of the black
line indicate insufficient velocity to stand, whereas points right
of the dashed line indicate a catastrophic fall forward. The
distance to the dashed line is used to measure stability, with
larger distances indicating higher stability.
Table V shows the median ROSv value across 5 trials for
each STS motion. In both cohorts, ROSv overall determines
that slow and quasi-static manoeuvres are more stable than the
fast and momentum-transfer manoeuvres, respectively, which
is congruent with the intuition developed in Section III-A.
Upon further examination, ROSv is most unreliable when
computing the stability of young subjects performing slow
and fast STS motions. Figure 5 shows the ROSv plot for a
young subject (a. ID 7) and an older subject (b. ID 11) that are
inaccurately characterised by ROSv. Hereafter, we continue to
highlight subjects ID 7 and 11 to compare the accuracy of each
method. These results illustrate the deficiencies of estimating
the stability of motion with a single feature and the inability
of the ROSv metric to characterise the specific perturbations
that lead to a fall.
Slow
Fast
Quasi-static
Momentum

1

0.8

0.6

0.4

0.2

0
-0.6

Slow
Fast
Quasi-static
Momentum

1

Normalized COM Velocity

m

Normalized COM Velocity

m2

IV. R ESULTS
In this section we compare an existing method of estimating
stability to the framework proposed above. All analysis was
performed on a system with an Intel Core i7-4930K 3.40GHz
processor with 12 cores and 32 GB RAM. The optimal
control problem was solved using MATLAB’s nonlinear solver
fmincon [55]. The optimization problem (D) is solved using
SPOTLESS [56] and MOSEK [57]. Code and figures may be
found at: https://www.w3id.org/people/vshia/jrsi.

0.8

0.6

0.4

0.2

-0.4

-0.2

0

0.2

0.4

0.6

0.8

Normalized COM Position

(a) Young Subject (ID 7)

1

0
-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Normalized COM Position

(b) Older subject (ID 11)

Fig. 5: ROSv plots of an individual’s STS motions. The stars denote individual
trials for each motion. The area between the black and dashed lines represent
initial positions and velocities where subjects have enough torque to stand.
The area to the right of the dashed line indicates falling forward, and the area
to the left of the black line indicates remaining seated.

B. Results for IPM
Using the method proposed in Section II, the BOS for each
individual’s motion is computed for the IPM. The optimal
5 Position

normalised to the subject’s foot length and velocity normalised to
pendulum length / second.

control was performed using 101 time steps with a polynomial
input of degree 6, and the optimization problem (D) was
solved with a degree 14 polynomial. The entire pipeline for a
single action required an average of 211 seconds to compute.
To compare the computed BOS for trajectories of different
time lengths, the BOS volume is normalised by the volume
of the domain with bounds. For example, 100% indicates that
the bounded domain is in the BOS and 0% indicates that the
BOS is empty. Table III shows the median of the normalised
volumes for the computed BOS of each subjects over all trials
of a specific manoeuvre.
All slower and quasi-static STS motions have larger basins
than faster and momentum-transfer STS motions, respectively,
indicating that subjects who use slower and more static motions are able to withstand more perturbations. The proposed
method correctly determines the STS strategies with higher
stability according to the intuition developed in Section III.
Group

Young

Older

ID
1
2
4
5
6
7
8
9
10
14
3
11
12
13
15

Untrained
Slow
35.9
42.9
38.5
38.1
41.0
31.2
37.3
35.4
36.1
37.7
38.8
40.3
35.5
51.2
45.9

Motion
Fast
23.0
17.2
16.8
18.4
9.7
16.5
21.9
18.5
10.7
16.7
26.1
21.0
19.1
23.4
16.8

Trained STS Strategy
Quasi-static
Momentum
34.5
11.9
46.9
15.9
30.4
13.9
38.1
19.3
45.1
8.7
37.5
15.0
37.5
17.2
38.7
14.3
38.9
8.8
35.3
11.6
42.2
23.4
45.5
14.5
34.6
21.1
53.0
9.9
40.6
12.0

TABLE III: Median volume of the BOS for each STS strategy calculated
using the IPM, as a percentage of the domain with bounds. The individuals
and strategies illustrated in Figure 5 are shown in bold.

The shape of the BOS in Figures 6 and 7 indicates the
perturbations the individuals in Figure 5 are able to withstand
under a specific control. Notice that the BOS for the quasistatic strategy is larger, indicating greater stability, at the
onset of the STS action. These results demonstrate that the
proposed method succeeds in cases that ROSv fails to properly
assess, and provides further information regarding sources of
instability throughout the STS action.
C. Results for DPM
To determine the effect of using a model that more accurately reflects the morphology of an individual, the BOS
for each individual’s motion is computed for the DPM. The
optimal control was performed using 101 time steps with a
polynomial input of degree 4, and the optimization problem
(D) was solved with a degree 8 polynomial. The entire
pipeline for a single action required an average of 4130
seconds to compute. Table IV describes the volume of the
BOS normalised by the volume of the bounded domain using
the the DPM for all subjects.
The BOS for subjects ID 7 and 11 are illustrated in Figures
8 and 9, respectively. Again, note for subject ID 7, the BOS

Group

ID
1
2
4
5
6
7
8
9
10
14
3
11
12
13
15

Young

Older

Untrained
Slow
16.9
21.1
13.1
12.9
17.8
18.4
13.5
27.6
14.8
23.0
36.5
28.1
8.5
18.0
24.1

Motion
Fast
11.4
11.0
10.2
8.2
16.3
12.2
7.5
17.4
9.3
14.7
30.6
20.2
5.9
15.4
12.0

Trained STS Strategy
Quasi-static
Momentum
20.8
11.1
30.6
17.0
14.6
10.8
15.2
8.4
22.8
13.2
19.5
12.6
13.8
7.0
26.3
19.3
17.0
11.3
24.6
11.2
32.7
29.8
27.6
22.1
9.8
6.7
26.3
14.8
19.6
11.3

TABLE IV: Median volume of the BOS for each STS strategy calculated
using the DPM, as a percentage of the domain with bounds. The individuals
and strategies illustrated in Fig. 5 are shown in bold.

for the slow STS is larger than the fast STS, and for subject ID
11, the BOS for quasi-static STS is larger than the momentumtransfer STS. Much like the IPM, the DPM also succeeds in
cases when ROSv fails. The BOS holds a rectangular shape at
certain times due to the θ1 and θ2 bounds and indicating that
the subject is “maximally” stable with respect to those states.
D. Summarizing Performance
As shown in Table V, the BOS computation method
presented above correctly identifies the STS strategies with
greater stability with a higher accuracy than the ROSv method.
Furthermore, by increasing model complexity from the IPM
to DPM, the volume of the computed BOS tends to decrease,
suggesting that higher order models may yield tighter BOSs
about the trajectory.
Method

Type

Fujimoto

ROSv
IPM
DPM

Proposed

Slow > Fast
Young
Older
6/10
4/5
10/10
5/5
10/10
5/5

Static > Momentum
Young
Older
9/10
4/5
10/10
5/5
10/10
5/5

TABLE V: The number of individuals for which each method correctly
identified the more stable standing strategy (Slow and Static). The median
value for ROSv and the BOS volume across 5 trials for each STS motion was
used for the comparison.

V. D ISCUSSION AND C ONCLUSION
This paper presents the first personalised computational
framework to model, identify, and analyse the stability of an
individual’s STS motion by using kinematic observations to
compute the BOS. Rather than reducing the STS motion to
a single feature, the entire trajectory is analysed via subjectspecific models and motion-specific trajectories to provide a
more informative metric of stability. Where ROSv method
fails, our proposed method successfully identifies slow and
quasi-static as more stable than fast and momentum-transfer
standing strategies. The shape of the computed BOS reveals
how stability changes throughout the STS motion, aiding in
the identification of unstable manoeuvres.

This framework provides a clinical tool to aid physical
therapists in identifying and reducing locomotor instability.
Using numerical tools to compute the BOS of locomotion
obviates the need to perform extensive perturbation experiments, making this approach applicable to injured or highrisk individuals. Double blind tests comparing the motion
of healthy, injured, and fall-prone individuals will help establish the diagnostic benefit of this automated computation
of stability. Because the shape of the BOS characterises the
perturbations most detrimental to an individual, this method
can direct the customization of physical therapy regimens
to improve stability. Furthermore, implementing our computational framework into longitudinal studies will improve
our understanding of the effect aging, injury, and clinical
intervention have on locomotor stability and quality of life.
Fortunately, the speed of this computation method makes it
feasible to quickly collect enough trials to experimentally
validate these clinical applications with adequate statistical
power.
There is an unavoidable tradeoff between computation
speed, accuracy and dimensionality. For systems with few
states (i.e. IPM), it is possible to simulate the volume of the
BOS via direct simulation and circumvent the optimization
problem defined in Section II-E. However, direct simulation
suffers from exponential scaling in the number of states.
For example, for the DPM, a 2-link pendulum, a sparse
simulation of the BOS consists of 175k randomly sampled
points requires over 8 hours to compute. To simulate the
BOS for a more representative human model such as a 3link pendulum or higher would take days or weeks, which is
not practically feasible for widespread deployment. However,
equally as important to the BOS is the identification of control
strategy used by humans for different actions.
This framework can be expanded to evaluate the stability
of a variety of human behaviours, thereby enabling the study
of human motion from a control-theoretic point of view.
Observations of athletes attempting to regain balance suggest
that swinging appendages can contribute to overall stability.
Although it has been demonstrated that the angular momentum
of swinging appendages can affect body rotation [58], the evaluation of control strategies exploiting swinging appendages
has not occurred.
The proposed method for computing the BOS enables
automation, widespread deployment, and customization of
motion analysis. This framework aids experimental design
and analysis of a variety of motions, enhancing the study
of individual differences in musculoskeletal architecture and
motor control strategies. The methods presented here can be
used to improve the identification of individuals at risk for
falling and to develop targeted therapy to increase stability,
thereby helping individuals maintain mobility and quality of
life.
VI. ACKNOWLEDGEMENTS
This work was supported by ONR MURI N00014-13-10341.

VII. C OMPETING I NTERESTS
The authors have no competing interests.
VIII. AUTHOR C ONTRIBUTIONS
VS formulated the framework and implementation, conducted human experiments and analysis, and drafted the
manuscript; TM assisted in data analysis and helped revise
the manuscript for important intellectual content; RB helped
conceive the study and draft the manuscript; RV aided in
the development of the framework, constructed the tools for
stability analysis, and helped draft the manuscript. All authors
gave final approval for publication.
R EFERENCES
[1] US CDC. Important Facts about Falls.
[2] Robertson MC, 2001 Effectiveness and economic evaluation of a nurse
delivered home exercise programme to prevent falls. 1: Randomised
controlled trial. BMJ 322:697–697. doi:10.1136/bmj.322.7288.697.
[3] Horak FB, 2006 Postural orientation and equilibrium: what do we need
to know about neural control of balance to prevent falls? Age and Ageing
35 Suppl 2:ii7–ii11. doi:10.1093/ageing/afl077.
[4] Pollock A, Durward B, Rowe P, Paul J, 2000 What is balance? Clinical
Rehabilitation 14:402–406. doi:10.1191/0269215500cr342oa.
[5] Campbell AJ, Borrie MJ, Spears GF, 1989 Risk factors for falls in
a community-based prospective study of people 70 years and older.
Journal of Gerontology 44:M112–7. doi:10.1093/geronj/44.4.M112.
[6] Berg K, 1989 Measuring balance in the elderly: preliminary development of an instrument. Physiotherapy Canada 41:304–311.
doi:10.3138/ptc.41.6.304.
[7] Lundin-Olsson L, Nyberg L, Gustafson Y, 1997 ”Stops walking when
talking” as a predictor of falls in elderly people. Lancet 349:617.
doi:10.1016/S0140-6736(97)24009-2.
[8] Podsiadlo D, Richardson S, 1991 The Timed ”Up & Go”: A Test
of Basic Functional Mobility for Frail Elderly Persons. Journal
of the American Geriatrics Society 39:142–148. doi:10.1111/j.15325415.1991.tb01616.x.
[9] Lundin-Olsson L, Nyberg L, Gustafson Y, 1998 Attention, Frailty,
and Falls: The Effect of a Manual Task on Basic Mobility. Journal
of the American Geriatrics Society 46:758–761. doi:10.1111/j.15325415.1998.tb03813.x.
[10] Pai YC, Patton J, 1997 Center of mass velocity-position predictions for balance control. Journal of Biomechanics 30:347–354.
doi:10.1016/S0021-9290(96)00165-0.
[11] Papa E, Cappozzo A, 1999 A telescopic inverted-pendulum model of the
musculo-skeletal system and its use for the analysis of the sit-to-stand
motor task. Journal of Biomechanics 32:1205–1212. doi:10.1016/S00219290(99)00103-7.
[12] Patton JL, Pai YC, Lee WA, 1999 Evaluation of a model that determines
the stability limits of dynamic balance. Gait & Posture 9:38–49.
doi:10.1016/S0966-6362(98)00037-X.
[13] Papa E, Cappozzo A, 2000 Sit-to-stand motor strategies investigated
in able-bodied young and elderly subjects. Journal of Biomechanics
33:1113–1122. doi:10.1016/S0021-9290(00)00046-4.
[14] Pai YC, Wening JD, Runtz EF, Iqbal K, Pavol MJ, 2003 Role of feedforward control of movement stability in reducing slip-related balance loss
and falls among older adults. Journal of Neurophysiology 90:755–62.
doi:10.1152/jn.01118.2002.
[15] Fujimoto M, Chou LS, 2013 Region of Stability Derived by Center
of Mass Acceleration Better Identifies Individuals with Difficulty in
Sit-to-Stand Movement. Annals of Biomedical Engineering 42:1–9.
doi:10.1007/s10439-013-0945-9.
[16] Bogle Thorbahn LD, Newton RA, 1996 Use of the Berg Balance Test
to Predict Falls in Elderly Persons. Physical Therapy 76:576–583.
[17] Steffen TM, Hacker TA, Mollinger L, 2002 Age- and Gender-Related
Test Performance in Community-Dwelling Elderly People: Six-Minute
Walk Test, Berg Balance Scale, Timed Up & Go Test, and Gait Speeds.
Physical Therapy 82:128–137.
[18] Muir SW, Berg K, Chesworth B, Speechley M, 2008 Use of the
Berg Balance Scale for predicting multiple falls in community-dwelling

[19]
[20]
[21]
[22]

[23]
[24]
[25]

[26]
[27]

[28]

[29]
[30]

[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
[40]
[41]
[42]

elderly people: a prospective study. Physical Therapy 88:449–59.
doi:10.2522/ptj.20070251.
Shumway-Cook A, Brauer S, Woollacott M, 2000 Predicting the Probability for Falls in Community-Dwelling Older Adults Using the Timed
Up & Go Test. Physical Therapy 80:896–903.
Prajna S, Jadbabaie A, 2004 Safety Verification of Hybrid Systems Using
Barrier Certificates. Hybrid Systems: Computation and Control pages
477–492.
Topcu U, Packard A, Seiler P, 2008 Local stability analysis using simulations and sum-of-squares programming. In Automatica, volume 44,
pages 2669–2675. IEEE. doi:10.1016/j.automatica.2008.03.010.
Mitchell IM, Bayen AM, Tomlin CJ, 2005 A time-dependent
Hamilton-Jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions on Automatic Control 50:947–957.
doi:10.1109/TAC.2005.851439.
Lyapunov AM, 1992 The general problem of the stability
of motion.
International Journal of Control 55:531–534.
doi:10.1080/00207179208934253.
Parrilo PA, 2000. Structured semidefinite programs and semialgebraic
geometry methods in robustness and optimization.
Maidens JN, Kaynama S, Mitchell IM, Oishi MMK, Dumont
GA, 2013 Lagrangian methods for approximating the viability
kernel in high-dimensional systems.
Automatica 49:2017–2029.
doi:doi:10.1016/j.automatica.2013.03.020.
Henrion D, Korda M, August D, 2013 Convex computation of the region
of attraction of polynomial control systems. IEEE Transactions on
Automatic Control 59:1–1. doi:10.1109/TAC.2013.2283095.
Majumdar A, Vasudevan R, Tobenkin MM, Tedrake R, 2014 Convex
optimization of nonlinear feedback controllers via occupation measures. The International Journal of Robotics Research 33:1209–1230.
doi:10.1177/0278364914528059.
Shia V, Vasudevan R, Bajcsy R, Tedrake R, 2014 Convex Computation of the Reachable Set for Controlled Polynomial Hybrid Systems. In IEEE Conference on Decision and Control. doi:DOI:
10.1109/CDC.2014.7039612.
Mohan S, Shia V, Vasudevan R, 2016 Convex Computation of the
Reachable Set for Hybrid Systems with Parametric Uncertainty page 25.
Koopman B, Grootenboer HJ, de Jongh HJ, 1995 An inverse dynamics model for the analysis, reconstruction and prediction of bipedal
walking. Journal of Biomechanics 28:1369–1376. doi:10.1016/00219290(94)00185-7.
Hargraves CR, Paris SW, 1987 Direct trajectory optimization using
nonlinear programming and collocation. Journal of Guidance, Control,
and Dynamics 10:338–342. doi:10.2514/3.20223.
Flash T, Hogan N, 1985 The coordination of arm movements: an experimentally confirmed mathematical model. The Journal of Neuroscience
5:1688–703.
Uno Y, Kawato M, Suzuki R, 1989 Formation and control of optimal
trajectory in human multijoint arm movement. Biological Cybernetics
61. doi:10.1007/BF00204593.
Kawato M, 1999 Internal models for motor control and trajectory planning. Current Opinion in Neurobiology 9:718–727. doi:10.1016/S09594388(99)00028-8.
Todorov E, Jordan MI, 2002 Optimal feedback control as a theory of motor coordination. Nature Neuroscience 5:1226–35. doi:10.1038/nn963.
Cusumano JP, Dingwell JB, 2013 Movement variability near goal equivalent manifolds: fluctuations, control, and model-based analysis. Human
movement science 32:899–923. doi:10.1016/j.humov.2013.07.019.
Liu D, Todorov E, 2007 Evidence for the flexible sensorimotor strategies
predicted by optimal feedback control. The Journal of Neuroscience
27:9354–68. doi:10.1523/JNEUROSCI.1110-06.2007.
Dingwell JB, John J, Cusumano JP, 2010 Do humans optimally exploit
redundancy to control step variability in walking? PLoS Computational
Biology 6:e1000856. doi:10.1371/journal.pcbi.1000856.
Dingwell JB, Smallwood RF, Cusumano JP, 2013 Trial-to-trial dynamics
and learning in a generalized, redundant reaching task. Journal of
Neurophysiology 109:225–37. doi:10.1152/jn.00951.2011.
Callier FM, Desoer CA, 1994 Linear System Theory. Springer New
York.
Rudin W, 1964 Principles of mathematical analysis, volume 3. McGrawHill New York.
Full RJ, Kubow T, Schmitt J, Holmes P, Koditschek D, 2002 Quantifying
dynamic stability and maneuverability in legged locomotion. Integrative
and comparative biology 42:149–57. doi:10.1093/icb/42.1.149.

[43] Spagna JC, Goldman DI, Lin PC, Koditschek DE, Full RJ, 2007
Distributed mechanical feedback in arthropods and robots simplifies
control of rapid running on challenging terrain. Bioinspiration &
biomimetics 2:9–18. doi:10.1088/1748-3182/2/1/002.
[44] Maki BE, McIlroy WE, 1997 The Role of Limb Movements in Maintaining Upright Stance: The ”Change-in-Support” Strategy. Physical
Therapy 77:488–507.
[45] Wand P, Prochazka A, Sontag KH, 1980 Neuromuscular responses to
gait perturbations in freely moving cats. Experimental Brain Research
38. doi:10.1007/BF00237937.
[46] Hof AL, Vermerris SM, Gjaltema WA, 2010 Balance responses to lateral
perturbations in human treadmill walking. The Journal of Experimental
Biology 213:2655–64. doi:10.1242/jeb.042572.
[47] Hughes MA, Weiner DK, Schenkman ML, Long RM, Studenski SA,
1994 Chair rise strategies in the elderly. Clinical Biomechanics 9:187–
92. doi:10.1016/0268-0033(94)90020-5.
[48] Aissaoui R, Dansereau J, 1999 Biomechanical analysis and modelling of sit to stand task: a literature review.
IEEE International Conference on Systems, Man, and Cybernetics 1:141–146.
doi:10.1109/ICSMC.1999.814072.
[49] Anan M, Ibara T, Kito N, Shinkoda K, 2012 The Clarification of the
Strategy during Sit-to-Stand Motion from the Standpoint of Mechanical
Energy Transfer. Journal of Physical Therapy Science 24:231–236.
doi:10.1589/jpts.24.231.
[50] AMTI. OPT464508.
[51] PhaseSpace Impulse X2. Technical report.
[52] PhaseSpace, 2010 Recap2 User’s Guide. Technical report, PhaseSpace.
[53] Wolfram Research. Double Pendulum.
[54] Pavol MJ, Owings TM, Grabiner MD, 2002 Body segment inertial
parameter estimation for the general population of older adults. Journal
of Biomechanics 35:707–712. doi:10.1016/S0021-9290(01)00250-0.
[55] The MathWorks, 2015 MATLAB: Optimization Toolbox User Guide.
Technical report, The MathWorks, Inc.
[56] SPOTLESS.
[57] APS Mosek. The MOSEK optimization software.
[58] Libby T, Moore TY, Chang-Siu E, Li D, Cohen DJ, Jusufi A, Full RJ,
2012 Tail-assisted pitch control in lizards, robots and dinosaurs. Nature
481:181–4. doi:10.1038/nature10710.

(a) Young subject - Slow

(b) Young subject - Fast

Fig. 6: The computed BOS for the IPM using the proposed method for untrained STS motions of subject ID 7. The top panels show a diagram of the STS
motion through time, with the non-gray colour depicting the IPM representation of the motion. The bottom panels show the computed BOS for the IPM in
the state space of the model. Thick coloured lines represent the observed trajectory (xobs ) of the STS motion. The gray region represents the computed BOS
with time slices in dark gray corresponding to the diagram above. Time is on the X-axis (scaled in each subplot), angle from foot to COM on the Y-axis,
and angular velocity on the Z-axis.

(a) Older subject - Quasi-static

(b) Older subject - Momentum Transfer

Fig. 7: The computed BOS for the IPM using the proposed method for trained STS motions of subject ID 11. The top panels show a diagram of the STS
motion through time, with the non-gray colour depicting the IPM representation of the motion. The bottom panels show the computed BOS for the IPM in
the state space of the model. Thick coloured lines represent the observed trajectory (xobs ) of the STS motion. The gray region represents the computed BOS
with time slices in dark gray corresponding to the diagram above. Time is on the X-axis (scaled in each subplot), angle from foot to COM on the Y-axis,
and angular velocity on the Z-axis.

(a) Young subject - Slow

(b) Young subject - Fast

Fig. 8: The computed BOS for the DPM for untrained STS motions of subject ID 7. The top panels illustrate the STS motion through time, with the coloured
skeleton depicting the DPM representation of the motion. The center and bottom panels represent the BOS for each joint of the DPM in the state space of
the model. Thick coloured lines in each subplot represent the observed trajectory (xobs ) of the STS motion. The gray region represents the computed BOS,
with time slices in dark gray corresponding to the position of the model above. Time is on the X-axis (scaled in each subplot), angle from the joint (1 or 2)
to the limb COM (m1 or m2 ) on the Y-axis, and corresponding joint angular velocity on the Z-axis.

(a) Older subject - Quasi-static

(b) Older subject - Momentum Transfer

Fig. 9: The computed BOS for the DPM for trained STS motions of subject ID 11. The top panels illustrate the STS motion through time, with the coloured
skeleton depicting the DPM representation of the motion. The center and bottom panels represent the BOS for each joint of the DPM in the state space of
the model. Thick coloured lines in each subplot represent the observed trajectory (xobs ) of the STS motion. The gray region represents the computed BOS,
with time slices in dark gray corresponding to the position of the model above. Time is on the X-axis (scaled in each subplot), angle from the joint (1 or 2)
to the limb COM (m1 or m2 ) on the Y-axis, and corresponding joint angular velocity on the Z-axis.

