main.tex

\documentclass{article}

% for images: png, pdf, etc
\usepackage{graphicx}

% for nice table formatting, i.e. /toprule, /midrule, etc
\usepackage{booktabs}

% to allow for \verb++ declarations in captions.
\usepackage{cprotect}

% for nice source code syntax highlighting, also provides Listing env
%\usepackage{minted}

% for nice units
\usepackage{siunitx}

% for less margin space
%\usepackage{fullpage}

% for nice todo notes
\usepackage{todonotes}

% macros
\newcommand{\mstar}{\mathbf{m}^*(\varphi)}
\newcommand{\K}{\mathbf{K}(\varphi)}

\author{Jason K. Moore and Antonie van den Bogert}

\title{Direct Control Identification in Human Gait}

\date{}

\begin{document}

\maketitle

\section*{TODO}

\listoftodos

\begin{abstract}
  \todo[inline]{Write one!}
\end{abstract}

\section*{Introduction}
%
Recent research and commercial activity indicate that gait-related powered
prosthetics will play an important role in both assisting humans with
disabilities and enhancing the abilities of the able bodied. These devices
include a variety of sensors and actuators than can be coupled to a control
system to provide assisted gait. There are currently commercial version of
powered trasnfemoral leg prostheses \cite{BiOM} and lower extremity exoskeleton
\cite{ReWalk} and many products still in the research phase
\cite{Ekso,Parker,Goldfarb,etc} that hope to increase the mobility of the
disabled. As promising as these devices are the available
lightweight lower extremity exoskeletons lack gait that resembles an
able-bodied human both in terms of stability and ``natural'' gait motions.

Due to the highly non-linear plant robust control is difficult to come by from
first principles for these devices. Although work based on \cite{Geyer2010} are
giving promising results in simulated musculoskeletal systems \cite{Wang2012,
Geitenbeek2013} these controllers still only produce marginally stable
simulations.

Even though constraining the controller to rudimentary models of the
musculoskeleteal system is fruitful, we are curious if more generic control
structures can reveal more robust controllers. To improve the gait of powered
prosthetics, our intent is to identify a simple, linear controller from a large
set of data collected from able-bodied subjects being perturbed by random
longitudinal forces.

We pose our problem as such: given a ``human-like'' plant, i.e. one with
similar degrees of freedom and inertial properties to a human, that is actuated
with simple joint torques can we identify a feedback control mechanism for this
system that causes it to walk stably and recover from perturbations in much the
same way a human does?

We hypothesize that typical gait measurements taken while a person is walking
under the influence of perturbations can be used to identify a generic and
potentially robust feedback controller for simulated or robotic system. Other
researchers are also pursing this path.

\cite{Sup2008} proposed a finite state proportional derivative feed back
controller for a powered transfemoral prosthesis. The control law that
generates the desired joint torques is based on three constants and a desired
angle set point for each of the four modes of a gait cycle. They found the
values of the controller by fitting their controller to normative joint angle,
angular rate, and torque estimates from a single gait cycle from an able bodied
walker. This method is unlikely to actually identify a true controller due
to the inherent errors expounded by \cite{Rouse2013} but may have merit as
\cite{Rouse2014} found quasi-stiffness to be equivalent to stiffness. The
authors ended up hand tuning their controller based on user feel leaving them
with a controller that was less than data driven.

\todo[inline]{Add in a description of Hugh Herr's BiOM controller}

Elliot Rouse's \cite{Rousse2014} work is most similar to ours. They identified
the parameters of a linear mass-spring-damper model of the lumped ankle's
passive and active impedance by perturbing at four points in the stance phase of the gait
cycle. He does not explicitly claim that their model is a controller, opting
to call it modulated impedance instead. Rouse et al. carefully isolate the
change in joint angle and torque due to the perturbation and derive their model
parameters to avoid the issues discussed in \cite{Rouse2013}. He shows a linear
relationship for stiffness with respect to 20-70\% phase.

More recently \cite{Wang2014} uses a simpler approach and is able to predict
foot placement with a simple linear model derived from kinematic data that
relates the mid-stance pelvic state to the following foot location. This is
questionable though, as there is no way to show that this is a control mechanism
and not simply a mechanical correlation. This analyses was performed without
external perturbations suggesting that the walker self-perturbs enough to glean
results.

We approach this ``control'' identification problem in a similar fashion as
\cite{Rouse2014} but by making use of a much richer data set that has the
potential to expose multi-input multi-ouput control mechanisms for the entire
body. We collect greater than 15,000 gait cycles from 11 subjects walking at
three nominal speeds on an instrumented treadmill where we measure the
kinematics of a full body marker set and a full set of ground reaction loads at
each foot while perturbing the subject by varying the belt speed with
Gaussian-like noise. We then identify a gait phase gain scheduled non-linear
controller with the so called direct identification method~\cite{Ljung1999}.

\section*{Methods}
%
\subsection*{Experiments}
%
The data includes a large number of gait cycles (> 19000) for 11 subjects
subjects walking at three nominal speeds,
\SIlist{0.8;1.2;1.6}{\meter\per\second}, on a treadmill with and without
pseudo-random longitudinal perturbations in the belt speed.
The data is a subset of the data reported in \cite{Moore2015} and we refer the
reader to that publication for the details about the experimental protocol and
data. The subject metadata parameters are given in Table~\ref{tab:subjects}.
%
\begin{table}
  \cprotect\caption{Information about the 11 study participants. The final
    three columns provide the trial numbers associated with each nominal
    treadmill speed. The measured mass is computed from the mean total vertical
    ground reaction force just after the calibration pose event. Generated by
    \verb|src/subject_table.py|.}
  \centering
  \input{tables/subjects.tex}
  \label{tab:subjects}
\end{table}

\subsection*{Data Preprocessing}
%
The ``raw'' data provided by \cite{Moore2015} consists of 100~\si{\hertz} time
series of marker and ground reaction loads in addition to event time
identifiers. We make use of an open source software package,
GaitAnalysisToolKit, to process the data. The processing follows these steps:
%
\begin{enumerate}
  \item Identify missing markers and the corresponding time instances.
  \item Replace missing markers with linearly interpolated values.
  \item Filter all signals with a forward/backward 2nd order low pass
    Butterworth filter with a cutoff frequency of 6~\si{\hertz}.
  \item Section the trials into unperturbed and perturbed sections based on the
    event times.
  \item Compute the planar inverse dynamics for the lower body: right and left
    ankle, knee, hip joint angles, joint angular rates, and joint torques
    \cite{Winter2009}.
  \item Identify the heelstrike times based on the vertical ground reaction
    force.
  \item Segment the time series into gait cycles based on the right foot's
    heelstrikes.
  \item Linearly interpolate the time series in each gait cycle to a specific
    number of data points so all gait cycles have the same number of samples.
  \item Remove any outlier gait cycles that were not properly identified due to
    abnormal ground reaction force measurements.
\end{enumerate}

After the initial processing above, the data is stored in three dimensional
arrays for each phase of the trial M x N x (q + p), where M is the number of
gait cycles, N is the number of phase instances in the gait cycle, and (q + p)
are the number of measurements. Figure~\ref{fig:angle-torque-comparison} shows
the mean and standard deviation of across the gait cycles for selected time
series for both unperturbed and perturbed walking. Furthermore,
Figure~\ref{fig:gait-cycle-stats-comparison} compares several gait statistics
of interest for unperturbed and perturbed walking. These two figures are
intended to show that there is variation in the subject's motion that is a
direct result of the longitudinal perturbations.
%
\begin{figure}
  \centering
  \includegraphics{figures/unperturbed-perturbed-comparison.pdf}
  \cprotect\caption{Right leg mean and $3\sigma$ (shaded) joint angles and
    torques from both unperturbed (red) and perturbed (blue) gait cycles
    from trial 20. We define the nominal configuration, i.e. all joint angles
    equal to zero, such that the vectors from the shoulder to the hip, the hip
    to the knee, the knee to the ankle, and the heel to the toe are all aligned.
    Produced by \verb|src/unperturbed_perturbed_comparison.py|.}
  \label{fig:angle-torque-comparison}
\end{figure}
%
\begin{figure}
  \centering
  \includegraphics{figures/unperturbed-perturbed-boxplot-comparison.pdf}
  % NOTE : Make sure to update these values (120, 519) if anything in the
  % script changes.
  \cprotect\caption{Box plots of the average belt speed, stride frequency,
    stride length, and stride width which compare 120 unperturbed (U: red) and
    519 perturbed (P: blue) gait cycles. The median is given with the box
    bounding the first and third quartiles and the whiskers bound the range of
    the data. Produced by \verb|src/unperturbed_perturbed_comparison.py|.}
  \label{fig:gait-cycle-stats-comparison}
\end{figure}

\subsection*{Controller}
\label{sec:controller}
%
We assume a black box~\footnote{black box in the system identification sense,
i.e. there are parameters in the model} control structure for the closed loop
system, Figure~\ref{fig:controller}. We first assume there is an unknown
non-linear plant, i.e. the open loop musculoskeletal system and treadmill, that
has an unknown time varying state. This plant can be perturbed by an external
exogenous input $w(t)$, in our case by varying belt speed, and is also driven
by joint torques $\mathbf{m}(t)$. An unknown sensor model generates feedback
signals, $\mathbf{s}(t)$, that are compared to a nominal gait phase dependent
reference trajectory $\mathbf{s}_0(\varphi)$ that describes the gait phase
varying state for cyclic gait when no feedback is necessary. The error in these
sensors are fed into the controller which generate corrective joint torques
that are added to the nominal gait phase dependent joint torques
$\mathbf{m}_0(\varphi)$ that correspond to the nominal trajectories,
$\mathbf{s}_0(\varphi)$. We then assume that when the actual state deviates
from the nominal state trajectory that the controller computes additive joint
torques based on the state trajectory error to compensate for variation in gait
away from the nominal.

The controller is designated as a gait phase scheduled gain matrix, $\K$, that
relates the error in the state trajectories to the additive joint torques. The
controller structure bounded by the box in Figure~\ref{fig:controller} can be
described by the following algebraic equation.
%
\begin{equation}
  \mathbf{m}(t) = \mathbf{m}_0(\varphi) + \mathbf{K}(\varphi) [\mathbf{s}_0(\varphi) - \mathbf{s}(t)]
  \label{eq:controller}
\end{equation}
%
\begin{figure}
  \centering
  \includegraphics{figures/control-system.pdf}
  \caption{The controller block diagram}
  \label{fig:controller}
\end{figure}

It is important to note that this control structure will effectively lump the
passive control gains with the active ones in addition to any other unmodeled
effects. Other authors refer to these gains as the impedance stiffness and
damping. This is desirable for our intentions, as we are only concerned with
mapping this identified controller to that of powered prosthetic device and it
is thus unnecessary to distinguish between these.

\subsection*{Direct Identification}
%
We employ a direct identification approach to find the optimal parameters for
the control structure described in Section~\ref{sec:controller}. We are aware
that direct identifcation has its drawbacks. It is well known that process
noise due to unmodeled affects in the identification model will bias the
results towards the inverse of the plant if the external perturbations are not
sufficiently large enough~\cite{Ljung1999,Kearney1990,Kooij2005}. Determining
if the perturbations are of large enough size is difficult because the process
noise is generally unknown.

\todo[inline]{I need to explain why we think our perturbations may be large
  enough in magnitude. The problem is that we get similiar identification
  results from both unpertrubed and perturbed data. THis can mean one of two
  things: (1) the person perturbs themselves enough and the identification is
  correct either way or (2) our perturbations are not enough, i.e. we'd get
  different results by perturbing more. The only thing that I can think of here
  is to show the necessary perturbation level for the quiet standing case and
  then say that we have similar magnitudes.}

We reformulate Equation~\ref{eq:controller} to make it amenable to linear
identifaction. The equation can be rewritten so that is linear in the gains and
a new term $\mathbf{m}^*$.
%
\begin{equation}
  \mathbf{m}(t) = \mathbf{m}^*(\varphi) - \mathbf{K}(\varphi) \mathbf{s}(t)
  \label{eq:linear-form}
\end{equation}
%
where
%
\begin{equation}
  \mathbf{m}^*(\varphi) = \mathbf{m}_0(\varphi) + \mathbf{K}(\varphi) \mathbf{s}_0(\varphi)
\end{equation}

Assuming that we can collect noisy measurements of $\mathbf{m}(t)$
and $\mathbf{s}(t)$, $\mathbf{m}^*(\varphi)$ and $\mathbf{K}(\varphi)$ are
identified simply by using linear least squares to regress the data. Given $m$
gait cycles with $n$ time steps in each gait cycle with $q$ controls and $p$
sensors there are $nq(p + 1)$ unknowns and $mnq$ equations so we need at least
$p + 1$ gait cycles to solve for the unknowns.

Details of coercing Equation~\ref{eq:linear-form} into the canonical form,
$\mathbf{Ax}=\mathbf{b}$ can be found in the supplementary IPython notebook
\verb|form_linear_system.ipynb|.

We compute the scheduled gains and $\mathbf{m}^*$ for both the unperturbed and
perturbed gait cycles for each trial. Three quarters of the data from each
series of gait cycles are used to compute the unknowns and the model is
validated against the remaining one quarter of the data. We compute the
percentatge of variance accounted for by the model with respect to the
validation data set to gauage how well the model predicts.

\begin{equation}
  \textrm{VAF} = 1 - \frac{||\mathbf{m}_p - \mathbf{m}_m||}
    {||\mathbf{m}_m - \bar{\mathbf{m}}_m||}
  \label{eq:vaf}
\end{equation}

\section*{Results}
%
The formulation above allows one to choose a variety of combinations of
potential sensor, $\mathbf{s}(t)$, and actuator, $\mathbf{m}(t)$, signals to
generate a controller. For this analyses, we assume we are constructing a
controller for a planar plant that only involves the hip, knee, and ankle
joints of each leg. We choose the angle and angular rate of each joint as
sensors and a torque at each joint are as the actuators.

The most general case for the controller, $\K$, has $q \times p$ unknown
entries per phase point but the formulation allows for different types of
control that isolate actuators from sensors. We explore four possibilities:
%
\begin{description}
  \item[Full Control] The acutation at a joint is affected by all of the
    available sensors.
  \item[Side Isolated Control] The actuation at a joint is only affected by the
    sensors on the same side of the body.
  \item[Joint Isolated Control] The actuation at a joint is only affected by
    the sensors at that joint.
  \item[No Feedback Control] The motion is only governed by $\mstar$.
    \todo{Isn't the ``no feedback'' equivalent to $\mstar$ being equivalent to
    the mean torques?}
\end{description}

%
% - 11 different subjects
% - qxp different mass normalized gains
% - n different points in the gait cycle
% - 3 different speeds
% - 2 genders
% - unpeturbed and perturbed
% - right and left sides
%
% - Are the mean subject gains at each speed significantly different than zero
%   at each phase point?
% - Are the stds wrt to subject in one speed large or small?
% - Do the mean subject gains at each speed vary wrt respect to gait phase
%   points?
% - Do the mean subject gains vary with respect to speed at each gait phase
%   point?
% - Are the mean subject gains significantly different between right and left
%   legs?

\subsection*{Example Result}
%
Once a controller is identified from a selection of gait cycles from a trial
the primary results of interest are the gait phase dependent gain values and
how well the model can predict the acutation of indepedent validation data. The
actual gain values may give some insight to the physiological nature of the
feedback mechanism. Here we present a typical result from a join isolated
controller structure. Figure~\ref{fig:example-gains} shows the estimates of the
mass normalized scheduled gains with respect to the percent gait cycle in each
leg as a function of the phase of the gait cycle.
%
\begin{figure}
  \centering
  \includegraphics{figures/example-identified-joint-isolated-gains.pdf}
  \caption{Example identified body mass normalized gait phase scheduled gains
    for the right (blue) and left (red) legs identified from a selection of
    gait cycles from a single trial. The top row shows the proportional gains,
    the middle row provides reference mean trajectories of the sensor
    measurements, and the final row shows the derivative gains. The columns
    correspond to the three joints for each leg. The shaded area bounding the
    gains gives the standard deviation of the parameters with respect to the
    fit variance. Gain values greater than zero indicate negative feedback and
    gains below zero indicate positive feedback.}
  \label{fig:example-gains}
\end{figure}

\todo[inline]{The mismatch from positive feedback to negative K values is
potetnially confusing, as Ton has pointed out before, but to me this is the
standard way to describe a feedback system in control literature.}

Figure~\ref{fig:example-fit} shows an example prediction of the joint torques
in the right and left legs by the identified control model based on independent
data from the same trial. \todo{This plot technically shows the VAF based on
the 20 samples per cycle comparison, yet I plot all of the measured values,
i.e. more like 100-120 samples per cycle.}
%
\begin{figure}
  \begin{center}
    \includegraphics{figures/example-identified-joint-isolated-fit.pdf}
    \caption{Predicted joint torques using the identified model shown in
      Figure~\ref{fig:example-gains} compared to those computed from inverse
      dynamics. The sensor and actuator estimations are taken from independent
      gait cycles from the same trial. The variance accounted for, VAF, is
      reported for all of the validation gait cycles period whereas only a select
      number of gait cycles are shown.}
    \label{fig:example-fit}
  \end{center}
\end{figure}

\subsection*{Inter subject mean}
%
We compute the body mass normalized gains, $\frac{\mathbf{K}(\varphi)}{m}$ for
each subject at each nominal speed and examine the mean and standard deviation
across subjects for the inter-subject differences.
Figure~\ref{fig:mean-gains-1-2} shows the mean and standard deviation of the
gains identified from a joint isolated control structure at the
1.2~\si{\meter\per\second} nominal speed.
%
\begin{figure}
  \begin{center}
    %\includegraphics{figures/example-mean-gains-1-2.pdf}
    \caption{Mean and standard deviation (shaded) body mass normalized gains
      across all subjects walking at a nomimnal speed of
      1.2~\si{\meter\per\second}.}
    \label{fig:mean-gains-1-2}
  \end{center}
\end{figure}


\subsection{Speed variation}
%
\section*{Discussion}
%
Figure~\ref{fig:example-gains} shows that there are patterns that are
significantly different than zero in portions of the gait phase for each gain.

\todo[inline]{How do I show that the gains are significantly different than
  zero?}

\subsection*{Gains are consistent across subjects}
%
Figure~\ref{fig:mean-gains1-2} indicates that the gains exhibit similar
patterns across subjects. There is substantial variability at some points in
the gait cycle. The variability is higher in 40\% to 50\%. \cite{Rouse2014}
found high intersubject variability in the ankle model parameters during the
stance phase.

\todo[inline]{Should I somehow check whether the subjects are significantly
  different from one another? i.e. check whether each subject is significantly
  different than the mean. Would it be useful to compute the coefficien to of
  variation (std/mean)? This could blow up for values around zero.}

\subsection*{No feedback in the swing phase}
%
During the swing phase of each leg the gains are not significantly different
than zero for the ankle and knee gains. This indicates that swing leg ankle and
knee may not be utilized in a feedback fashion to react from the perturbations
orginating in the stance leg.

\todo[inline]{It would be nice to plot the mean toe off time on the gain plots
  to delineate the stance and swing phase. In addition I could add the gait man
  to the plots.}

\subsection*{Gains exhibit anatomical symmetry}
%
The gait data was split into gait cycles with respect to the right leg's
heelstrikes. The gains identified from the right delineated gait cycles exhibit
anatomical symmetry with respect to the right and left legs.

\todo[inline]{It is probably possible to compute whethter the right is
significantly different than the left.}

\subsection*{Directly Identified Gait Phase Scheduled Feedback Gains Exhibit
  Positive and Negative Feedback}
%
The derivative gains exhibit a combination of negative and positive feedback.
For example, in the middle of the stance phase the ankle derivative gain
indicates a large positive feedback. This implies that if the foot deviates
from the nominal angular rate, the controller will push it

If the ankle angle deviates from the nominal angle the torque is such that you
will be pushed back toward the nominal position, yet at the same time if the
angular rate deviates from the nominal

\todo[inline]{This is really hard for me to think about what happens.}

\subsection*{Random variations predict random gains}
%
If the variation in the sensor and actuator measurements is purely random, the
identified gains should be random. To verify this we follow this procedure:

\begin{itemize}
  \item Select the necessary marker coordinates and ground reaction load
    measurments for a single gait cycle from normal walking.
  \item Compute the inverse dynamics from sensor and actuator measurements.
  \item Fit a Fourier series to each time series so that the cycle can be
    perfectly cyclic.
  \item Form a time series of 500 gait cycles based on the fitted data.
  \item Add normal Gaussian noise to each sensor and actuator measurement.
  \item Identify a joint isolated controller from this artificially created
    data.
\end{itemize}

The identified gains show no apparent correlation. \todo{Add a gain plot here
showing this}

\subsection*{Correlations due to inverse dynamics bias gains}
%
The measurement noise present in the marker coordinates is propogated into both
the sensor values (joint angles and angular rates) and the actuator values
(joint torques) through the inverse kinematic and inverse dynamic computations,
thus the sensors and actuator measurement estimates are correlated simply by
this. We exposed this correlation by the following procedure:

\begin{itemize}
  \item Select the necessary marker coordinates and ground reaction load
    measurments for a single gait cycle from normal walking.
  \item Fit a Fourier series to each time series so that the cycle can be
    perfectly cyclic.
  \item Form a time series of 500 gait cycles based on the fitted data.
  \item Add normal Gaussian noise to each measurement.
  \item Compute the inverse dynamics from the noisy data to produce the sensor
    and actuator measurements.
  \item Identify a joint isolated controller from this artificially created
    data.
\end{itemize}

The correlation based entirely on the random noise is shown in the example gain
plot in Figure~\ref{fig:inverse-dynamics=correlation-gains}.
%
\begin{figure}
  \begin{center}
    \includegraphics{figures/example-inverse-dynamics-correlation-gains.pdf}
    \caption{Gait phase percent scheduled gains for right (blue) and left (red) legs.}
    \label{fig:inverse-dynamics=correlation-gains}
  \end{center}
\end{figure}

The proportional gains show a clear bias that is strikingly similar to a
proportionally scaled vertical ground reaction force.

\subsection*{$\mathbf{m}^*$ accounts for most of the variation}
%
Ideally the variation in the gait cycles is sufficient to identify the
acutation due to feedback encapsulated by the gait matrix. If the feedback
contribution is assumed to be zero, then $\mathbf{m}^*$ can be identified
alone. If the VAF decreases significantly with $\mathbf{K}{\phi}=0$ then one
can assume that the feedback terms are not just contributing to overfitting the
data.
%
\begin{table}
  \cprotect\caption{Generated by \verb|src/table_no_control_vaf_comparison.py|.}
  \centering
  \input{tables/no-control-vaf-table.tex}
  \label{tab:no-control-vaf}
\end{table}

Table~\ref{tab:no-control-vaf} shows that $\mathbf{m}^*$ alone is sufficient to
explain a large portion of the variance in the data. Adding the feedback terms
pushes the VAF closer to 100\%.

\subsection*{Similar gains are predicted from both unperturbed and perturbed}
%
This means that either that there is enough feedback variation in normal
walking to identify the controller or that are not identifying anything that is
significant.

\section*{Conclusion}
%
We are able to identify a simple linear controller that exhibits larger gains
in the stance phase than in the swing phase. Additionally, similar gain
patterns in the right and left legs are observed that use both positive and
negative feedback. The controller is capable of predicting the measured joint
torques with greater than 65\% VAF in all joints. Results and conclusions from
a larger sample of subjects and conditions will be presented at the conference.

\section*{Acknowledgments}

This research was funded by the Ohio's Wright Center for Sensor Systems
Engineering.

\bibliographystyle{plain}
\bibliography{references}

\end{document}