5"!”

,‘1

‘1

 

r ..

Jam

37

 

-.
:52“;
E“

 

THESIS

illMINillHIWlilllllllllilllxHiWWI

3 1293 01050 41

 

LIBRARY

Michigan State
University

 

 

 

This is to certify that the

dissertation entitled

9051.14.01? Feta/[aft ELS/IMa/IbV)

P/&Ceo€bvr€ «7/? film" (Jed/C
M45Ccl/oékc/cila/ Yél'C n4

!‘ '1 fMdC
via 6 presentedb§l Kan/an 251%! 4V)

WIS/AI Y OHfiH ER

has been accepted towards fulfillment
of the requirements for

13H . 0 degree in [7/ €C/) 6(1); :5

@m 5- was

Major professor

Date [2-7, 76

MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

PLACE IN RETURN BOX to romovo this checkout from your record.
TO AVOID FINES return on or odor. duo duo.

DATE DUE DATE DUE DATE DUE

      

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU lo An Affirmative Adlai/Equal Oppoflunlty lnotltmon
W

ulna-9.1

 

 

 

 

POSITION FEEDBACK ESTIMATION
PROCEDURE IN A LARGE SCALE
MUSCULOSKELETAL SYSTEM
VIA EXTENDED LINEARIZATION

By

Yasin Yousef Dhaher

A DISSERTATION

Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Materials Science and Mechanics

1 996

ABSTRACT

Position Feedback Estimation
Procedure in a Large Scale
Musculoskeletal System
Via Extended Linearization

By

Yasin Y. Dhaher

Synthetic studies of human cybernetics have been the focus of intensive studies in
the past decade. Considered a branch of motor control, sensory feedback, modulated by
sensory information, responds to changes in system states. The nature of the feedback
adaptability is considered to be fundamental to the functional coordination in human
motion activity.

This study proposes an approach to compute the position feedback in an arbitrary
human activity. Pivotal to the method is the idea that any human activity can be
characterized by an input/output mathematical model where the system inputs and
outputs are the muscular activity and the motion nominal trajectories, respectively. This
model includes both the human dynamics and control elements, and therefore is a closed
loop system. The parameters needed for the closed loop model are obtained using the
predictive error parametric identification method. To describe human dynamics, an
anthropomorphic musculoskeletal model is proposed that includes the same number of

inputs and outputs measured in the identification of the closed loop system. Since the

anthropomorphic model is free of any feedback components, it is defined as the open
loop system. The open loop model is nonlinear, thus an extended linearization method
was used. The linearization procedure reduced the nonlinear model to a family of
linearized systems parametrized by the same set of measured operating points used in the
estimation of the closed loop system. A key step in the procedure is to represent both the
estimated closed loop system and the proposed open loop model in the same state space
representation. The block observer form was utilized as the standard state space
realization of the closed and open loop systems.

Once the open and closed loop systems are available, the feedback gain estimation
is obtained by finding the feedback gain matrix that when combined with the proposed
open loop model resembles the estimated closed loop structure. An explicit solution of
the feedback gain matrix and a detailed discussion of the methods to validate the solution
are included. The proposed method is demonstrated for the movement of squatting.
Quantitative results of the positional feedback gains are presented. Results of agonistl
antagonist position feedback synergy in the muscles included in the model are also
shown. A brief discussion on the use of the estimation method proposed is given in the

areas of predictive orthopaedic surgeries and neurological deficiencies.

COpyright by
Yasin Y. Dhaher
1996

 

Acknowledgments

I take this opportunity to show my appreciation to the following people whose contri-

butions made the completion of the thesis possible

To my major professor, Dr. Robert Soutas-Little, for his continued support and encour-

agement.

To Dr. Hasan Khalil, for sharing his expertise in the development of the control theory

part of the present work.
To my committee members: Dr. Gary Cloud, Dr. Joseph Vorro, and Dr. Robert Hub-

bard.

To Stefanie Bonin, Patricia Soutas Little, LeAnn Slicer and all my peers (past and

present) at BEL.

iii

 

 

Table of Contents

List of Tables

vi

 

List of Figures

 

Chapter One
Introduction

vii

 

1.1 Why Study The Feedback Problem

 

1.2 Background, Modeling, Identification and Control
1.3 Overview of the Thesis

 

 

Chapter Two
Theoretical Development

 

2.1 Identification of the Closed Loop System

 

2.2 Modeling of the Open Loop System

 

2.2.1 Musculotendon Moment Arm Computation

 

2.2.2 Muscle Modeling and Architecture

 

2.2.3 Musculoskeletal Mathematical Model

 

2.3 Estimation of the Position Feedback Structure

 

Chapter Three
An Illustrative Example

 

3.1 Anthropomorphic Musculoskeletal Model, ‘Open Loop’

3.2 Closed Loop Model

 

 

3.2.1 Experimental Setup

 

 

3.2.2 Experimental Results
3.2.3 Autoregressive Closed Loop Estimation

 

3.3 Position Feedback Solution

 

15

l 9
20
29
29

33
41

45

54
54
68
68
72
79
83

3.3.1 Validity of the Solution of the Position Feedback Problem .................. 83

iv

3.3.2 Interpretation of Some of the Position Feedback Elements ................. 87
Chapter Four

 

 

 

Conclusions 106
Appendices

Appendix A - 110
Appendix B 113
Appendix C 115

 

References 120

List of Tables

 

Table l. Musculotendon static properties of the muscles used in the
analysis. 62
Table 2. The body-segmental parameters used for the skeletal model. ............. 69

Figure 1.1

Figure 2.1.1

Figure 2.1.2A

List of Figures

Conceptual representation of the human global control as proposed

 

 

 

by Bernstein, 1947. 13
Schematic illustration of the identification procedure defined in the

present study, where ts is the sampling period. 25
A typical electromyographic signal. 28

Figure 2.1.2B A comparison between the output of the bilinear model of Zajac,

Figure 2.2.1

Figure 2.2.2

Figure 2.2.3

Figure 2.3.1

Figure 3.1.1

Figure 3.1.2

Figure 3.1.3

Figure 3.1.4
Figure 3.1.5

Figure 3.1.6

1989 and the signal processing method used in the present study.

 

 

28
Schematic sketch of an m-articular muscle used for muscle moment
arm and length calculations. 30
Schematic sketch of the muscle architectural model developed in
Gordon et al., 1986. 34

 

Simulation display of the muscle nonlinear contraction dynamics
of the Rectus Femoris and hamstring muscles. The rate of tension
defined in equation (2.2.10) is computed as shown for both muscles.

 

 

 

 

 

40
Schematic illustration of the output feedback estimation procedure
developed in this study 46
The sagittal plane linkage model including a set of muscle groups
representing the squat activity under study. 55
The sagittal plane linkage model showing the segmental coordinate
systems. _ - 59
Maximum isometric ankle plantar flexion moments, with knee in full
extension. 64
Maximum isometric ankle dorsiflexion moment. 64

 

Maximum isometric knee extension moments, with hip in full
extension. 66
Maximum isometric hip extension moments, with knee in full

vii

Figure 3.1.7
Figure 3.1.8
Figure 3.2.1
Figure 3.2.2
Figure 3.2.3
Figure 3.2.4
Figure 3.2.5

Figure 3.2.6

Figure 3.2.7

Figure 3.2.3

Figure 3.2.9

Figure 3.3.]

Figure 3.3.2

Figure 3.3.3

Figure 3.3.4

Figure 3.3.5

Figure 3.3.6

Figure 3.3.7

extension.

 

Maximum isometric hip extension moments, with knee in full
extension.

 

Maximum isometric flexor moment of the rectus femoris muscle
about the hip joint, with knee in full extension
An illustrative sketch of the targeting scheme used in the present
study on the shank and the location of the laboratory coordinate
system in the calibration space.

 

The joint coordinate system of the shank body segment ......................

Segmental angles.

 

Muscle activity plots of hamstrings and gluteus maximus for
three squat cycles

 

Muscle activity plots of rectus femoris and the vasti muscles for
three squat cycles.

 

Muscle activity plots of the gastrocenimus muscle for three
squat cycles.

 

Muscle activity plots of the soleus and the tibialis anterior muscles

for three squat cycles muscle activity plots of the gastrocenimus
muscle for three squat cycles.

 

Measured segmental angles (solid lines) and simulated

results according to the MIMO ARX model (dashed lines). .............

Stick figures of the four degrees of freedom rigid body segment

in squat activity. The solid lines are the measured motion data, the

dashed lines are simulated data.

 

Schematic representation of the validation procedure.

 

Relative percentage error plot for the computation of the shank
angle.

 

Relative percentage error plot for the computation of the trunk
angle.

 

The stick figure and the associated joint motions for the one squat

cycle

 

Feedback gain plot from the ankle, knee, and hip motions to the
tibialis anterior muscle.

 

Feedback gain plot from the knee and the hip motions to the
soleus muscle.

 

Feedback gain plot from the ankle and the hip motions to the

viii

67

67

7 1
74
76

78

78

79

82

86

86

89

vastus muscle- 91
Figure 3.3.8 Ankle angle feedback gain plot to the soleus muscle for a squat

 

 

 

 

 

activity. 93
Figure 3.3.9 Ankle angle feedback gain plot to the tibialis anterior muscle for a

squat activity. 94
Figure 3.3.10 Knee angle feedback gain plot to the rectus femoris muscle for a

squat activity. 96
Figure 3.3.11 Knee angle feedback gain plot to the rectus femoris muscle for a

squat activity. 97
Figure 3.3.12 Knee angle feedback gain plot to the hamstrings muscle for a

squat activity- 99
Figure 3.3.13 Hip angle feedback gain plot to the hamstrings muscle for a

squat activity. 100

 

Figure 3.3.14 Position feedback gain associated with the ankle angle cross plot
of the soleus and tibialis anterior muscles for a squat activity.
102

 

Figure 3.3.15 Position feedback gain associated with the knee angle cross plot
of the hamstring and rectus femoris muscles for a squat activity.
103

 

Figure 3.3.16 Position feedback gain associated with the hip angle cross plot
of the gluteus maximus and the hamstring muscles for a squat
activity. 105

 

ix

CHAPTER ONE

INTRODUCTION

1.1 Why Study the Feedback Problem?

The complexity of the human body requires an advanced control system. This
control system continuously manipulates and processes sensory information relayed by
local and global sensory components. In this manner, sensory feedback is considered as
the principal component in movement organization and control in humans. Stability,
precision, and hierarchical control of movements are considered to be the major
products of the sensory system. Since general human activity involves large motions,
the nature of a feedback system is self-tuning. Responding to positions, velocities, and
forces, the feedback system reschedules its sensory information to obtain the three
products as defined above. The relationships between physical quantities (position,
velocities, forces. etc.) and the electrical signals produced by the sensory system are
known as feedback gains.

Classical approaches used to compute these feedback gains, such as the state-
feedback analyses, fail to address the adaptability of such gain values. State-feedback
methods, for example, are useful in the study of human activities that exhibit linear
dynamics such as postural control. It is also important to note that the state-feedback
analysis depends on the assumed values of the closed loop eigenvalues. The sign and the

magnitudes of these eigenvalues are always chosen to maintain a stable closed loop

system. Unfortunately, classical methods give no physical justification for the choice of
these values.

Since complex motions, for example, rising from a chair, are mathematically repre-
sented in terms of nonlinear dynamics, a self-tuning feedback control structure is required
to accommodate the large changes in system configurations. Thus, the current constant
feedback studies fail to represent the muscular synergies of the agonist and antagonist
actuators controlling the activity.

The development of a method to compute the feedback gains that deals with the non-
linear dynamics of the open loop system as well as the actual dynamic behavior of the
closed 100p system, provides significant challenges in the field of human cybemetics.
This type of analysis has not been published in the literature. The proposed procedure to
compute the feedback gains depends on the actual open and closed loop systems. Funda-
mental to the method is the construction of input / output mathematical models for both
the closed and the open loop systems. The open loop structure is obtained using theories
in musculoskeletal modeling. The nonlinear form of the musculoskeletal model is con-
verted to a family of parameterizeparameterized set of linear systems using extended lin- ‘
earization approach. System identification theory is used to construct an input / output
closed loop mathematical model from data measured on a live subject.

The computation of feedback gains are important to areas of the closed loop control
design of the functional electrical stimulation systems (Stanic and kaoczya, 1974; Kralj
et al., 1980; Crago and Chizeck, 1986; Khogh and Zajac, 1989, MD. Applications of
the feedback analyses are essential to understanding the sensory components’ contribution

to the stability, precision and control of the human motion (Brain, 1989; Iqbal and

Hemami, 1993; Kuo, 1995). The application of computing the feedback gains to the field
of predictive orthopaedic surgeries would prove to be beneficial in future research.

1.2 Background, Modeling, Identification and Control

In the following discussion, the literature on concepts used in the development of the
theory of this thesis will be presented. To follow the development of the theoretical
model, the discussion starts with an overview on the theory of system identification and
its application to the field of biodynamics. Next, a detailed presentation of basic research
and ideas concerning modeling of musculoskeletal structures are given. Lastly, the theo-
ries and research on human control system are discussed.

The dynamics of a system can be described by a mathematical model that relates its
inputs and outputs. The inputs and outputs are a set of time sequences that represent the
external variables acting on the system, as well as the observable responses of that sys-
tem. In general, the system can be investigated either theoretically or experimentally.
The theoretical analysis of a process is based on the use of a set of physical balance equa-
tions and the understanding of the phenomenological laws. Experimental analysis treats
the system as a black-box with a mathematical structure that consists of a set of elements
to be estimated through the use of any of the available identification routines (Ljung,
1987; Zhang et al, 1990). The sequential input/output relation that describes the process
dynamics is defined as the transfer function represented in terms of a backward shift oper-
ator (q) as used in control theory (Strejc, 1981). In the multi-inputs multi-outputs
(MIMO) case, the transfer function will not be a single function, but a matrix that in gen-

eral takes the form

111(9) hm”) hlm(qfl
H”) = h21(q) h22(g) h2m(q)

 

 

th1(g)hp2(q) hpm(q)
where hij(g) are in general rational functions of q. The input/output relation can then

be written as follows

Y(q) = H(q) ' U(q)
Where Y(q) and U(q) are the output and input vectors, respectively. The problem of

finding the coefficients of the functions hij(q) from measured input and output data is com-

monly known as transfer function identification (Soderstrom et al., 1989). The identifica-
tion of the transfer matrix can be done using either parametric or non-parametric
identification procedures (Jenkins et al., 1969; Soderstorm, 1984; Ljung, 1986). In the
non-parametric identification method, the identification is carried out without employing
a set of parameter vectors in search for the best description, hence the name non-para-
metric (Jenkins et al., 1969). Unlike non-parametric identification, parametric identifica-
tion methods do not depend on test signals to identify a model. A selected model is
chosen and is characterized by a set of parameters. The best set of parameters to mini-
rnize a selected criterion is obtained by an estimation process. The prediction error
method is one of the most commonly used parametric identification procedures. The
search process in the prediction error method is based on the least square method of linear
regression where the parameter vectors are found such that the sum of the squares of the
output prediction errors are small (Lawson and Hanson, 1974; Ljung, 1986). The choice

of the parameters is highly dependent on the model chosen. Models are usually selected

from a wide range of already established models (Soderstorm, 1984). The Autoregressive
(ARX) and the Autoregressive Moving Average (ARXMA) models are among the com-
monly used system structures (Ljung, 1986).

Treating the electromyographic myosignals as system inputs and the joint or segmen-
tal angles as system outputs, the human musculoskeletal system transfer function can be
modeled and identified. Zhang et al., (1990) considered the modeling and identification
of human knee joint dynamics. Eight myosignals representing eight muscles crossing the
knee joint were considered as the system inputs. Data for the six degrees of freedom of
the knee joint were selected as the system outputs. Several methods were then tested to
carry out the identification procedure of a proposed linear, (ARX), and nonlinear model.
The results obtained using the prediction error method for both the linear and nonlinear
models were comparable. Other identification methods, like the instrumental variable
method and the non-parametric identification method, were also used to study human
postural global control (Ishida and Miyazaki, 1987; Johansson and Magnusson, 1991).

Modeling and identification of human structures are important to understand the syn-
thesis and control of the musculoskeletal system. Identification, for example, can be
useful in areas such as the restoration of functional tasks to paralyzed muscle—limb sys-
tems by means of the functional electrical stimulation (FES) (Franklen et al., 1995; Fran-
klen et al., 1995). Franklen at al. (1995) examined a system which consisted of the
quadriceps electrically stimulated using surface electrodes. A parametric model was then
derived to adequately represent the active component of the muscle that would be used to
enhance the control performance of the electrically stimulated paralyzed limbs.

Theoretical mechanical models of the segmental human musculoskeletal body appear

frequently in the literature (Nubar and Contini, 1961; Beckett and Chang, 1968; Kane
and Scher, 1970; Krogrnan and Johnston, 1970; Von Gierke, 1971; Chao and Rim, 1973;
Hatze, 1973; Seirge and Arvikar, 1973; Hatze, 1976; and others). The first analytical
model that defines human dynamics was introduced by Fisher in 1906. Fisher developed
a set of equations of motion of an n—links model using Lagrangian dynamics.

Before the introduction of methods and the development of an analytical dynamic
model, some of the difficulties that are inherent characteristics of the musculoskeletal
system will be explored. Owing to the complexity of the human musculoskeletal system,
(144 joints, including the minor joints, and approximately 750 muscles) assumptions
should be made to reduce the number of joints and muscles that are to be included in the
model. The inclusion of all the six degrees of freedom of each joint in a model results in
about 850 degrees of freedom to represent the whole system dynamics. Clearly, the inclu-
sion of about 850 degrees of freedom combined with about 750 muscles results in an
overly complicated system. Hence, simplifications have to be introduced in the process
of modeling the skeletal system. The question is: How many muscles and joints and
- degrees of freedom of each joint should be included in the model? Decisions are usually
made depending on the type of activity considered in the study. For example, in the case
of squatting, the dominant motion of all the joints is in the sagittal plane; and, in the
absence of pathology, a one degree of freedom model of the included joints may be ade-
quate (Hemami, 1978; Hatze, 1981; Huoston and Passarello, 1982).

The other major component of the musculoskeletal model is the muscle that acts as the
actuator for force generation to drive the skeletal system. Two main issues are important

when considering muscles in a model; the muscle architecture; the moment arm about the

joint that the muscle spans. Although Fisher’s work is considered to be the foundation of
studies of human mechanics, his work was insufficient to explain the generalized forces in
the Lagrange’s formulation. Since the moment arm is fundamentally one of the compo.
nents in the computation of the generalized forces, Hatze (1965) was the first to introduce
a transformation based method through which the moment arm of a muscle about a joint
can be represented in terms of that joint’s angle, provided that the muscle origin and
insertion data are known. The method developed by Hatze was modified by a more recent
work presented by Gordon et al. (1986). In both the early work of Hatze and the most
recent work by Gordon et al., the anthropometric data of the muscle origin and insertion
are of most importance. Thus, many cadaver based studies, focusing on obtaining mus-
cle origin and insertion points were done (Alexander and Vernon, 1975; Brand and
Crowninshield, 1982; White et al., 1989; Seirg and Arviker, 1989; Komistek et al.,
1994; and others).

Inclusion of a muscle model in the analysis of human dynamics is significant since dif-
ferent muscle models affect the prediction of muscle forces, hence joint moments and
generalized forces. Mathematical modeling of muscles has been .the subject of research on
both the molecular (see Needham (1971) for a review) and macroscopic (Carlson, 1957;
Green, 1969; Crowe, 1970; Pierrynowski and Morrison, 1985) levels. Macroscopically,
muscles with tendons are divided into two major parts, the muscle tendon and the muscle
belly. The tendon, a collagenous tissue structure, has been shown to exhibit a nonlinear
viscoelastic behavior (Fung, 1967; Haut and Little, 1972; Soong and Huang, 1973; Jhon-
son et al., 1992). The muscle belly consists of many individual muscle fibers held

together by connective tissue (fascia). This fascia is called the epimysin. Another type

of fascia is the perimysin. The perimysin is a fascia that penetrates the muscle belly sepa-
rating muscle fibers into groups called fasciculi. Each fasciculus consists of a number of
muscle fibers wrapped with another connective tissue known as the endomysin. Collec-
tively, the epirnysin, perimysin, and endomysin act as passive elements that are structur-
ally located parallel to the muscle fibers, hence the name parallel elements. D.K. Hill -
(1968) showed that most of the tension observed when stretching a resting muscle comes
from the parallel elements. Based on quick release tests, Soong and Huang (1973)
showed that a constitutive equation of the muscle parallel element takes a nonlinear expo-
nential form. The presence of muscle series elements was theorized to exist based on
experiments formed by Hill A. V. (1938; 1950; 1953); Mlkie (1956); Sonnenblick (1964);
and others. Unlike the parallel element, the muscle series element exhibits elastic behav-
ior that is activation dependent. Functionally, the series element is responsible for the
transmission of tension to the end points of the muscle fiber when stimulation takes place
(Hatze, 1975). To model a constitutive equation of the series elastic element, Pier-
rynowski and Morrison, 1985, developed a model using experimental data from earlier
research done by Bahler (1967); Bahler and Fales (1968); and Close (1972). More details
on the muscle series element are found in the comprehensive study of Ehema and Huijing
(1990).

The only force-generating activation dependent muscle element is the muscle contrac-
tile element. Structurally, the contractile element represents a set of cascaded building
blocks known as sarcomeres. Sarcomere shortening develops due to the sliding (cross-
bridging) of its major components, actin and myosin filaments. The contractile force pro-

duced by a muscle fiber is equal to the sum of all cross-bridges in one half-sarcomere of a

fiber. The functional contribution of the contractile element to the muscle machine is its
force-velocity relation. A.V. Hill (1938) was one of the first to propose an experimentally
based model relating the contractile element force with its shortening velocity. A rectan-
gular hyperbolic relation was developed and is known as Hill’s model. Other models
were also developed to characterize the actin/myosin bridge force generation and contrac-
tion. All of the studies discussed so far suffered from the inability to explain the stretching
behavior of the muscle model (Pringle, 1960). To account for the lengthening of muscle
fibers, Sugi (1972) modified Hill’s model based on his experimental observation of the
muscle fiber force generation in the presence of stretching. He observed that, with the
increasing stretching velocities, the force at first rises above the isometric value but then
levels off. Most recently, Hatze (1990) provided a comprehensive discussion on the phys-
ical and experimental complications in modeling the contractile element force-velocity
relation.

One may wonder if all the elements of the muscle should be included in the model and
how complex the mathematical model of every element should be. To answer these ques-
tions the model chosen should pass the following test: can the model predict the muscle
force generated in both muscle’s lengthening and shortening phases. It should also be
noted that increased complexity of the muscle model will subsequently increase the com-
putational time for analysis of the complete human system dynamic synthesis. Thus, the
dynamics synthesis computational time is also a factor that affects the choice of the mus-
cle model. Winters and Stark (1987) provided a good summary of the effect of the muscle
model selection on the amount of information gained or lost in the process of human

dynamic studies. A non-dimensional generalized model that may be used in computer

10

simulation was first introduced by Zajac et al. (1986). This model is dependent on the
muscle static properties, tendon slack length, muscle maximum isometric force, and
other factors. Details on the development of this model are given in the review by Zajac
(1989).

Once skeletal and muscle models are chosen, the derivation of the musculoskeletal
dynamic equations is commonly done by Lagrangian dynamics (Goldstein, 1959; Maros
and Orlanda, 1971; Leu and Hemati, 1986; Vukobrtatovic et al., 1990), Euler-Newton’s
(Chace and Bayazitolque, 1971; Wittenburg, 1977), Hamiltonian (Vane and Sitchin,
1970; Hagedom, 1979), and others. Lieh et al. (1990) provided an informative review of
the existing multi-body dynamics formulation techniques and their application to biome-
chanical systems.

The dynamics of human movement is controlled by inputs from sensory systems, or
triggered by sensory signals or some internal desire to produce a movement. The two sys-
tems interact in a way so that rhythmic movement is monitored and continually com-
manded for reinforcement and accuracy (Herman et al. 1976). Relevant to the present
work, the sensory feedback mechanism is two parts, a local feedback component defined
as a spinal segmental loop, and a global feedback component that is modulated through a
higher level of the central nervous system (CNS) as defined in Brooks (1989). The main
local and global feedback structures in the musculoskeletal system are the joint sensors
(Boyd and Roberts, 1953; Mlliams, 1981; Johansson and Magnusson, 1991), proprio-
ceptive sensors in the joint capsules, tendons, skin and muscles (Houk er al., 1970; Agar-
wal et al., 1970; Boyd, 1980; Johansson and Magnusson, 1990), vestibular (Nashner,

1971; Horak et al., 1990; Barker, 1991), and visual (Reichardt and Paggio, 1979; Rei-

11

chardt, 1980; Barker, 1991) systems. Both joint and proprioceptive sensors consists of
receptors in joints and muscle tissues. These sensors detect the relative position, motion
and forces at the joints and in the muscle tendons (Magnus, 1926; Grigg and Greenspan,
1977; Hasan, 1983). Attempts were made to model and describe the general transfer
function of these sensors (Grigg and Greenspan, 1977; Diener and Dichgans, 1988).
However, a complete model of the joint and proprioceptive sensors that can be used for
the analysis of human dynamic structure does not yet exist (\Vrlliams, 1981; Hemami,
1985). On the other hand, both the visual and vestibular systems were studied intensively
and reported in the literature. The vestibular system provides both the positional and the
dynamic information needed to stabilize the system (Young, 1970; N ashner, 1973; Mag-
nusson, 1986). The vestibular system responds to dynamic stimulation, (linear and angu-
lar acceleration), through its sensory organs in the semicircular canals. However, the
positional information that the vestibular system provides is due to the interconnections
and integration of information provided by the visual system. For example, the vestibular
system provides postural stability and updates the system with spatial awareness by pro-
viding an eye tracking system such that a visual contact is maintained during head move-
ment. Functional models were developed to simulate the vestibular systems and their
contribution to postural stability (Young, 1970; Nashner, 1973). Awareness of position is
mainly attributed to the visual system. The visual feedback provides the system with
information concerning the relative position between a body and external references (Leh-
man and Stark, 1983).

A large number of analytical studies have been conducted to model and study the

effects of sensory mechanisms on human postural control and stability. Nashner (1971,

12

1972) developed a model of the postural balance in a rigid body of one degree of freedom
about the ankle joint. The focus of this study was the investigation of postural control
strategies and their relation to the proposed model. Later models, similar to those hypoth-
esized by Camana et al. (1977) and Golliday and Hemarni (1976) were proposed to repre-
sent both the velocity and the position feedback sensory modalities in a biped postural
model. He et al. (1991) developed a complex model where they included the feedback
pathways that represents the upper motor neurons to account for both the vestibular and
visual systems. They also included other sensory modalities intrinsic and extrinsic to the
muscle. Utilizing optimal control theory, they investigated the relationship of the selec-
tion of the input and states weighting matrices on the effect of the sensory components
proposed in the postural stability problem. Iqbal and Hemarni (1993) explored a concep-
tual feedback model of the proprioceptive sensory systems and studied the effect of the
model on the sway stability in a four link model.

An alternative approach to the micro-modeling of the sensory components is the
black-box approach. Analytically, the black-box method lumps all sensory feedback
mechanisms of the same nature in one element. Further decomposition of the estimated or
computed elements is then made to account for the different modalities (Hemarni and Gol-
liday, 1977; Hemarni and Jaswa, 1978; Brain, 1989; Kuo, 1995).

The control of a human activity is the most challenging issue to be resolved. About
fourty years after the first modeling attempt done by Fisher, Bernstein in 1947 investi-
gated the global control of human dynamics. He suggested that the overall control of the
musculoskeletal system consists of six components: 1. actuator (muscle), 2. sensor

(receptor) that senses information from both the skeletal and the muscular systems, 3. a

l3

programing device (CNS) which describes the necessary value of controlling parameters
to the system, 4. discriminator which gives the difference between the desired and actual
control parameter, 5. a coder system that decodes the signals from the discriminator and
sends them to the regulator (controller) through a feedback loop, and finally, 6. the con-
troller that controls the actuator function. A conceptual illustration of Bemstein’s inter-

pretation of the human motion control is shown in Figure 1.1.

 

 

Controller

 

 

1 '“'mg Discriminator

l W

Input
w —~ M... m

Human
Activity

Figure 1.1 Conceptual representation of the human global control
as proposed by Bernstein, 1947.

 

 

 

 

 

 

 

 

However, the criterion by which the controller (regulator) makes up the control sig-
nals was not clearly identified. Theoretical studies were then developed to give more
understanding of the form by which the control procedure defined by Bernstein is carried
out. Optimal control (Chow and Jacobson, 1971; Hatze, 1976; Khang and Zajac, 1989;
Pandy et al., 1990; Kuo, 1995), parametric optimization (Pandy et al., 1992; 1995),

constant state feedback (Zheng and Hemami, 1984 ; Brian, 1989), force feedback (Whit-

l4

ney, 1977; Raibert and Graig. 1981), nonlinear feedback (Hemami and Camana, 1976),
and Lyapunov stability (Hemami and Cvetkovic, 1977; Takeyaki and Arimoto, 1981;
Iqbal and Hemami, 1993) theories were used to define the control procedure in the human
structural control.

The basic question of the optimal control theory is to find a control input that mini-
mizes a predefined performance criterion with and without terminal conditions (Owens,
1981). In general, the solution of the optimal control problem leads to the solution of a
two point boundary value nonlinear problem defined by a set of equations known as the
Reccati equations (Owens, 1981). The size and the degree of the nonlinearity of the Rec-
cati equations depends on both the number of degrees of freedom included and the number
of input signals introduced in the system. The coefficient matrices in the Reccati equa-
tions are defined in term of input and output weighting matrices. Ho (1991) proposed
different forms of these matrices to account for different control strategies in the postural
sway problem. He assumed that for the linear postural system, the muscle excitations are
independent, hence the input weighting matrix was chosen to be diagonal. On the other
hand, the state weighting matrix was determined by the control strategies simulated. It is ‘
interesting to note that both the form and values of the matrices are of no physical signifi-
cance.

An alternative method of solving the optimal control problem is to convert the problem
to a parametric optimization problem (Sirisena and Tan, 1974; Goh and Teo, 1988;
Nagurka, 1990). Nagurka and You (1990) showed that by expanding the generalized
coordinates of the system in terms of Fourier series, the optimal control problem reduces

to finding the coefficients of the series expansion that would minimize a given objective

15

function. At every iteration, the inverse-dynamic equations are solved for the set of the
control inputs, hence the name inverse parametric optirrrization method. The iterations
continue until a minimum value of the objective function is reached (Gill and Murray,
1981). The major disadvantage of this method is that it can not cope with optimal control
problems that are Bang-Bang where the control input takes its limiting values and
switches between these values at given time intervals (Nagurka and Yen, 1990; Pandy et
al., 1992). Unlike the indirect method, in the direct parametric optimization procedure
described by Goh and Teo (1988), the system direct dynamic equations are solved for
every iteration seeking a minimal value of the objective function. The procedure starts by
parametrizing the input controls in terms of nodal points and then the system dynamic
equations are solved at discreet points in time. At every point, the objective function is
evaluated and compared with its previous value. Pandy er al. (1992, 1995) and Tashman
et al. (1995) applied the direct optimization method to a large-scale musculoskeletal
model. The direct method involves a large CPU time usage in the processes of reaching an
optimal solution. A detailed discussion of the CPU implications of the direct optimization
method is given by Ziegler et al. (1992). In both theories, optimal control and parametric
optimization, a set of control inputs are the main outcome of the analysis and they are not
constructed to compute any form of feedback structure that may represent the internal
feedback of the system. Studies involving the computation of the feedback structure
found in the literature, including the ones mentioned above, only deal with an activity

that exhibits a linear dynamic model, as in the case of the postural control problem.

1.3 Overview of the Thesis

Approaches to compute the position feedback gains in fields of human dynamics and

16

robotics involve either the use of optimal control theories or the state feedback analysis.
In general, the optimal control problem reduces to the solution of a set of nonlinear first
order differential equations representing the necessary optimal conditions. Unfortunately,
closed form solutions of these equations are hard to obtain due to the mixed boundary con-
ditions on the states and co-states equations. Therefore, numerical solutions are found by
applying forward integration to the state equations and backward integration to the co-
state equations. Numerical solutions are time-consuming, thus they are not robust
enough in cases where the feedback structure is self-tuning. It is also well known that the
complexity of the optimal control problem increases with highly nonlinear system dynam-
ics, as in the case of a large scale human activity. The state feedback analysis, on the
other hand, is generally based on the assumptions made over the closed loop eigenvalues
or characteristic equation. Since robots are man-made structures, the closed loop eigen—
values are usually chosen to meet certain performance requirements. However, in human
activities, less is known about these eigenvalues, beyond the fact that they are stable.
This study describes an alternate approach for computing the position feedback gains
in a large scale human activity. Pivotal to the method is the idea that any human activity
can be characterized by an input/output mathematical model where the system inputs and
outputs are the muscular activity and the nominal motion trajectories, respectively. Since
the input/output data are measured from a human subject, the model represents both
human dynamics and control, hence is defined as the closed loop system. The identifica-
tion of such a model from input/output data is obtained by using parametric identification
methods. The estimation procedure and issues related to the input and output data are dis-

cussed in Section 2.1.

17

To describe the human dynamic model, an anthropomorphic skeletal model that
includes the same number of inputs and outputs used in the identification of the closed
loop model is proposed. Implications, such as the number of degrees of freedom that
best define the activity under study, discussed in the previous section should also be
addressed in the modeling of the skeletal structure. A detailed development of the skeletal
model is found in Sections 2.2.1 and 2.2.3. This model is generally nonlinear, thus a lin-
earization procedure is invoked. Utilizing the extended linearization method proposed in
the literature, see Rugh (1984), the nonlinear dynamic equations of the skeletal model
are represented in terms of a family of linear systems parameterizeparameterized in terms
of the structure input/output operating points. The dynamic system’s nominal output data
are those measured and used in the identification of the closed loop system. However,
since the inputs of the skeletal dynamic system are the muscle forces, the system’s nonri-
nal input data are not available for measurements. Analytical methods are used instead.
Section 2.2.2 discusses an analytical method to compute the nominal muscle force based
on a muscle mechanics model given in the literature concerning the relationship of the
muscle force trajectories and the corresponding muscle activation curves. Both the com-
puted and measured input/output data together with the dynamic model are combined and
manipulated to form the family of linearized equations given in Section 2.2.3. Since the
anthropomorphic model is free of any feedback components, it is defined as the open loop
system.

The next step in the analysis is to estimate the position feedback gain matrix from the
available closed and open loop models. A key point in the procedure is to represent both

the open and close loop structures in the same state space model. The block observer state

18

space model is used as the standard structure for the analysis of this thesis. The computa-
tion of the feedback gain matrix is based on the following question: what would the feed-
back matrix be such that when combined with the open loop model it resembles the
closed loop structure? The procedure and condition of good estimation are given in Sec-
tion 2.3.

Finally, in Chapter 3, the utility of the proposed method is demonstrated by applying
it to the squat activity. Also in that chapter are qualitative results of the feedback coeffi-

cients and results of agonist/antagonist synergies of the muscles included in the model are

shown.

CHAPTER TWO

THEORETICAL DEVELOPNIENT

This chapter presents a quantitative method to predict the position feedback for a given
human activity. The method is based on two major structures, the open and closed loop
systems. As stated in the previous chapter, the proposed procedure starts by noting that
for given inputs (muscular activity) and outputs (ioint’s or segmental angles), a mathemat-
ical model that defines the input/output relation can be identified. Since this model is a
representation of the particular activity of a live subject, it will include both the subject
dynamics and control thus defined as the closed loop structure. In the present work, a
parametric estimation procedure was utilized to estimate the closed loop model. The esti-
mation procedure used is presented in Section 2.1.

The next step in the analysis is to propose an anthropomorphic skeletal model that best
represents the activity under study and also includes the same number of inputs and out-
puts used in the estimation of the closed loop system. Since the proposed model is free of
any form of feedback, it is defined as the open loop model. A detailed discussion of the
development and modeling of the open loop model is given in Section 2.2. The resulting
open loop model (skeletal linkage model) is generally a nonlinear model, so an extended
linearization method is invoked. The linearization procedure described in Section 2.2.3
uses the framework introduced by Rugh (1984) and is based on the idea that the nonlinear

open loop system can be replaced by a family of linearized systems parameterizeparame-

19

20

terized by a set of joint (segment) motions and nominal muscle forces data (operating
points). The nominal motion data are commonly measured using either photogrammetric
techniques or goniometers. On the other hand, a noninvasive nominal muscle force mea-
suring technique does not exist. Thus, alternate analytical methods are used to compute
the nominal muscle forces during a given activity. Among these procedures is the method
of static optimization (Seireg and Arivkar, 1973; Penrod et al., 1974; Crowninsheld and
Brand, 1981). Crowninshied and Brand used the static optimization theory by utilizing a
nonlinear performance criterion based on muscle endurance. One of the disadvantages of
using the optimization method is the large computational time needed to compute these
forces. However, Crowninshied and Brand found that the muscle nominal force com-
puted was in agreement with the linear envelope of its electrical activity (\Vrnter et al.,
1980). Consistent results were also found recently by Yamaguchi (1990). Thus, in the
present work, the physiological model developed by Gordon et al., (1986) is used to eval-
uate muscle nominal forces computed from measured muscle activities. The proposed
method is dependent on the understanding of a muscle phenomenological model used,
thus a detailed discussion of modeling a muscle architecture is given in Section 2.2.2.
Finally, Section 2.3 presents the position feedback estimation procedure used when both
the open and closed loop systems are mathematically defined.

2.1 Identification of the Closed Loop System

The terms “closed” and “open” loop systems will be used frequently in this work. The
closed loop system is defined as the actual dynamics of the human neuromusculoskeletal
system characterized by a set of experimentally obtained myosignals and segmental

angles. The dynamic activity observed is a product of not only the pure rigid body dynam-

21

ics but also includes the neurological system feedback and control. The open loop system
is defined as a mathematical model that consists of a set of articulated rigid bodies con-
nected by joints.

The goal of the present study is to construct a position feedback structure that, when
combined with a proposed open loop mathematical model, resembles an experimentally
developed closed loop structure. The computation of such a feedback structure depends
on the existence of a high quality closed loop model of the dynamics and control. The
main purpose of system identification theory is to develop a mathematical model that
defines the physical system from input/output experimental data. Various identification
algorithms, such as least squares, maximum likelihood, instrumental variable method,
cross correlation, and stochastic approximation have been applied successfully to the
parametric identification problem. Experimentally, the input/output data are available

only at discrete points in time to, t1, t2, . . . These instances of time can be arranged as
integral multiple of some basic unit ts, say 0 ts, 1 ts, 2 ts. . . In which case ts is often

known as the sampling period and the instants of time are defined in terms of a time
parameter k, where k takes the values 0, 1, 2, 3, . . . The sequential input/output (I/O)

relations that describe the process dynamics are known as difference equations. A typical

single-input single-output (8180) 11th order difference equation takes the form (Strejc,

1981)

an'Y(k‘")+an—r 'Y(k-(1-n))+...+y(k) _-.

b0 - u(k)-i-bl - u(k- 1) + +bn - u(k—n) +e(k) (2.1.1)

where y(k) is the output (generalized coordinates), u(k) is the input (activations), and

22

e(k) is the error associated with either the model or the experimental data. A delay oper-

ator q, is defined such that

(q-l)n-u(k) = u(k-n) n = 1,2,... (2.1.2)
The delay operator 9 has the same role in discrete systems that the Laplace transform
has in continous systems. For example, the g—transform of Equation (2.1.1) leads to an

algebraic equation in terms of the delay operator g. The resulting relation is known as the

transfer function between the input and the output. Invoking the delay transform of Equa-
tions (2.1.2) on Equation (2.1.1) gives

-(n-1)

arr—“flaw..-” °y(q)+...+y(q) =

b0 - u(q) +1:l - q“ -u(q) + +12, . u(q) +e(q) (2.1.3)
01'

""'”+...+1)-y(q) =

(an ° q-n+an—l ° 9
(b0 ”)1 - g" + +bn-g'")u(q)+e(g) (2.1.4)

In the multi-input multi-output (MIMO) case, Equation (2.1.4) can be written as
Am - Y(q) = 8m - th) +E(q) (2,1,5)

where Y (g) e R' is the output vector, U ( g) e R" is the input vector, A(q) and
B(q) are polynomial matrices of q with the appropriate dimensions. In Equation (2.1.4),
A(q) and B( 9) were (1x1) matrices where

A(g)=an.q‘"+a,,_,-q""'”+...+1 (2.1.6)

23

-n

B(g)=bo+b,.g“+...+b -q (2.1.7)

It

For the multi-variable model of Equation (2. 1.5), A( q) and B( q) generally take the fol-

lowing form

A(g)=An.q"’+An_l-q""'”+...+1 (2.1.8)

7'

—n

B(q)=Bo+Bl-q'l+...+Bnoq (2.1.9)

where A land B,- are constant matrices of appropriate dimensions and I r is an r x r identity

matrix. It will become apparent in the following sections that, for the purpose of the
present study, the order of the difference Equation (2.1.5), It takes the value 2. The
model developed and described above is commonly called a Discrete Autoregressive
(DARX) (MIMO) (I/O)-model (Isermann, 1989).

The identification process is the process of finding the appropriate coefficient matrices

of Equations (2.1.8) and (2.1.9), A], A2, ..., Bo, Bl, ..., that will satisfy a certain crite-

rion. The prediction error method is used in this study, as it is able to deal with the para-
metric identification of MIMO systems. The predictive error method is one of the least
square based algorithms that has been used successfully in the parametric identification of
problems that involve the parametric estimation in a large scale musculoskeletal systems
(Zhang et al, . 1990; Franklen et al., 1995). A detailed study of the method appears in the
work of Ljung (1986).

Rewriting Equation (2.1.5) in its equivalent difference equations the following is

obtained

An"’<k-n>+A...i-I/(Ic-n+1)+...+y(k) =

24

Bo - U(k)+Bl - U(k— 1) + +8" - U(k—n)+E(k) (2.1.10)

Define R as

R = [111,112 ...,13,,,19l ...] (2.1.11)
Then Equation (2. 1.10) can be rewritten in the following form

Y(k) = R'-‘P(k)+E(k) (2.1.12)

where

t
‘1’“) = [—Y'Uc— 1) ,—Y'(k—2) U‘(k) ,U‘(k— 1) ...] (2-1-13)

and

Y(k—i) = [y1(k-i) J2(k—i) yp(k-i)]t i = 1,2...n (2.1.14)

The form given in Equation (2.1.12) is known as the regression equation. The next
step in the estimation procedure is to solve Equation (2.1. 12) for the parameter matrix R.
A complete study of the available methods developed to accurately compute R from Equa-
tion (2.1.12) is presented in Ljung (1989). Based on the predictive error method, the esti-
mation of R involves an iterative checking scheme where the difference between the actual
output and predicted output is examined for a minimum value. The study of the conver-
gence and precision of the method is available in the system identification literature (Iser-
mann, 1989; Ljung, 1986; Soderstrom and Stocia, 1989; etc).

The block diagram shown in Figure 2.1. illustrates the general concept of the identifi-
cation procedure proposed here for the closed loop system. The availability of the exper-
imental data at discrete points in time is represented on Figure 2.1.1. In this study, the

closed loop inputs are considered to be the muscle activations and the outputs are the sys-

25

tem measured generalized coordinates. A comprehensive discussion on the nature of the
system inputs is given in this section, while the system outputs are discussed in detail in

Chapter 3.

 

 

 

 

 

 

 

Muscles neralized
Actrvatrons oordrn ates

{ Inputs I A human activity {outputs }

D ‘Neuromusculoskeletal’ * >

System
tS
Identification
procedure

 

 

 

Figure 2.1.1 Schematic illustration of the identification procedure
defined in the present study, where ts is the sampling period.

For the purpose of the current study, the inputs were the myosignals. The muscle
force generation process can be generally divided into two major mechanisms, muscle
activation and muscle contraction. The later will be discussed in detail in the following
sections. Muscle activation is defined by the electrochemical process that takes place
upon the arrival of action potentials (AP) at the neuromuscular junction at the terminal

arbor. Ebashi and Endo (1968) stated that the activation dynamics of a muscle are best
described by the concentration of the Ca” ions in the intera-filamentary space where the

activation state is defined by the relative amount of Ca” ions bound to the troponin
(Ebashi and Endo, 1968; Hatze, 1980; Pierrynowski and Morrison, 1985; Shipping and

Zahalak, 1988). If the maximum number of potential interactive sites on the actin fila-

26

ments are exposed by the action of the calcium, then the activation is said to be maximum
and takes a value of one. The muscle initial activation is defined as the minimum poten-
tial interactive sites that are present while the muscle is at rest. Hatze (1977) proposed a
second order lumped system that represents the hypothesis first proposed by Ebashi and
Endo (1968) and later was found to be consistent with the microscopic model developed
by Shipping and Zahalak (1988). Their model considered the action potential signal as the
system input (myoelectrical signal). Levine and Zajac (1984) developed a first order bilin-
ear contraction dynamic model as an activation model simpler than that developed by

Hatze (1980). The general form of the bilinear model as described in Zajac (1989) is

%am = {$113 + (1 - 13) - u(m] - am +§ - u(t) (2.1-15)

where u(t) is the input neural electromyosignal, a(t) is the muscle activation, 1 is the
activation process time constant, and B is a constant parameter defined as the ratio of the
muscle activation time constant over the muscle relaxation time constant, 0 < B < 1.
Burker et al. (1973) experimentally found that the activation time constant is the same for
the three types of muscle fibers, slow twitch, fast twitch and fast twitch fatigable fibers.
They estimated the activation time constant to be approximately 0.003 seconds. Typical
values of B ranges from 1/3 to 1/2 depending on the type of the muscle fiber (Close,
1972; Altringham and Jhunson, 1982). As a demonstration, Equation (2.1.15) was simu-
lated with an input that represented a typical electromyographic (EMG) signal shown in
Figure 2.1.2A. As shown in Figure 2.1.2B, the output of the bilinear model of Equation
(2.1.15) basically rectifies and modulates the raw EMG signal. Hence, in this study an
equivalent signal processing approach, that obtains the same outcome as that of the bilin-

ear model, was used. The experimental EMG signal was processed by first filtering and

27

then rectifying and modulating. A band-pass Butterworth filter of 5H2 and 200Hz cut-off
frequencies was used to filter the EMG. Then a full wave rectifier accompanied with a
low-pass Butterworth filter, 8Hz cut-off frequency, was used to obtain the corresponding
muscle activation curve. The result of the procedure is then considered to be the system
input used in the identification procedure presented earlier in this section. This signal pro-
cessing procedure is shown to be equivalent to the output of the model of Equation
(2.1.15), see Figure 2.1.28, and was recommended by Shiavi et al., 1985; Zajac, 1989;
\Vrnter, 1990; Zhang, 1990 and others as an alternative to solving and simulating Equa-

tion (2.1.15).

15

o
h

Normalized EMG
In 0

I
‘

L
in

h

‘

Normalized Activation
O
in

28

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_Raw Signal i
)— ..................... ' ooooooooooooooooooooo i ........ 1. .......... é ....... .Ip
I. I‘ II N l E 1' .14 ‘ II .‘I 1 ‘ L

1 l l l I i l L

O 01 02 03 04 05 01 07 OD

 

 

 

 

 

 

 

 

——Bilineer Model Output :3
... ....... t .......... .3 ........... 3.... ....... .; ..... _
. . 3 if; !
---- Srgnal Processmg Output , 5 3, 5
t E El E
I E l 3' 3
I . I 3
l : , E
+- ........................ e.- . .e ' e no '...........: ........ .0... ‘ 1 o u d
r I r :1 I - 3 ‘
I n 2 l 2 2 .
s r 1 ,' r r ‘ EH ’ , s s s)
‘ H a I E L J: \r E "
V 1" 1 i 1 i 4' 1
O 01 02 03 DA 05 0‘ 0] 03

Time (mseconds)

Figure 2.1.2A. A typical electromyographic signal.
Figure 2.1.2B. A comparison between the output of
the bilinear model of Zajac, 1989 and
the signal processing method used in
the present study.

29

2.2 Modeling of The Open Loop System

The open loop system, as defined in the previous section, consists of a system of
articulated rigid bodies. The choice of the number of muscle forces and the degrees of
freedom of the rigid bodies or the joints, is dictated by the number of inputs and outputs
observed in the identification procedure presented in Section 2.1. For example, if the
knee joint is modeled with six degrees of freedom, then all the six quantities should be
available for measurement. Muscles in musculoskeletal structures can be modeled as
either moment generators or sub-systems (Hatze, 1980). Including muscles as sub-sys-
tems will only increase the number of the states and equivalently increase the order of the
system. However, the conceptual implementation of the proposed position feedback esti-
mation procedure will not change. In this study, muscles are modeled as moment genera-
tors. In the next sub-section, a brief representation of the muscle moment arm
computation is presented. Then a complete discussion on muscle modeling and architec-
ture is presented in Section 2.2.2. Finally, the rigid body skeletal system mechanical
modeling is presented in Section 2.2.3.

2.2.1. Musculotendon Moment Arm Computation

To illustrate the computation of a muscle moment arm, consider the sketch shown in

-> ->
Figure 2.2.1. In the sketch shown, r 0 and r , represent the vectors to the origin and inser-

tion points with respect to their body segment coordinate systems, respectively. Yamag-

—> ->
guchi et al. (1990), provided tabulated values of r0 and r, for a large number of human

muscles from experimental studies on cadavers. Generally, a muscle may span more than

one joint and overlap more than one body segment (the hamstrings, for example). There-

 

 

 

 

in, h:

where n is

men! the 1

Mus

Imus

30

fore, Figure 2.2.1 illustrates a general muscle spanning (m) joints and (11) body segments,

where n is equal to (m+1). In Figure 2.2.1 the coordinates (xi, 1’], Zj), j =1, 2, . . n, rep-

resent the n segmental coordinates of the n body segments the muscle spans

    

Y n
f \ . Xn
r \ .
l \ ~
mJ'
Origin (0) "
9
r0 Z
n
Muscle Belly /
Tendon
, 9
/x R0
(I) Y,
Insertion
ii
I I XI
Zr

Figure 2.2.1 Schematic sketch of an m-articular muscle used for muscle moment arm and
length calculations.

It can be easily shown that the position vector of the origin point with respect to the

31

insertion point coordinate system (X1, Y], Z I) is given in its matrix form as

m m j
R0 = uni-0+ 2{HTJ._1 .pj} (2.2.1)

where T]- is the transformation matrix of (X14 1, Yj+ 1, ZJ-H) coordinates to (Xj, Yj, Zj)

P.

coordinates and represents the rotational degrees of freedom of j'h joint, J - J is the

1’
position vector of the (X14 1, 13+ 1, 21-”) coordinate system with respect to the (Xj, Yj, Zj)

coordinate system defined in the “I“ Y}- Zj) coordinates and represented in its matrix
form, and To is an identity transformation matrix. In general, P} is not function of time
unless the joint is said to have translational degrees of freedom. In that case, Pi is taken
to be a variable function of time that represents some fixed value and defines the structural
distance between the consecutive joints, and the translational J]- ‘h joint degrees of free-
dom. Once the position vector of the origin point R0 is computed, the insertion to origin
position vector is given as

——-> -) —>
RO/l = Ro-rl (2.2.2)

In cases where the muscle may have more than one insertion or origin points, (for
example the tibialis anterior muscle), effective origin or insertion points are introduced,
(see Gordon et al., 1986; Hoy et al., 1990) and the formulation above holds true.

The moment arm of the muscle can be computed from Equation (2.2.2) above using
the fact that the muscle moment arm is always perpendicular to the line of action of the

muscle force as shown in Figure 2.2.1. Thus, the moment arm of the muscle about for

example, the rum joint Jm , is given as

32

\
I

m
—->
m]. = (H Tj’rOJXRo/l/
i=1

—->

R0,, (2.2.3)

 

 

It should be obvious from Equation (2.2.3) that the muscle moment arm is a function

—->
of system degrees of freedom through both R o / I and T j‘s. This dependency adds another

complexity to the musculoskeletal system through the introduction of coupling in the

dynamic equations as will be shown in later sections. It can be easily shown that the

._)
moment of the muscle, 1," , about the J", joint shown in Figure 2.2.1 is given as

\
I

 

 

 

m R
o/I
-—) [H Ti ° r a] p R—_)
T] = p j = 1 1 x 0/1
m f \—.W_—J
moment arm muscle force
vector vector (2.2.4)

where p is the muscle force computed from a proposed muscle model (see the next sec-
tion for details).

Finally, both the muscle moment and moment arms computed using Equations (2.2.3
and 2.2.4) may also be used in cases where the joint exhibits translational degrees of free-
dom. In these cases, the origin point position vector i: will consists of two parts. The
first represents the anthropomorphic origin point location, obtained from the literature,
and the second component characterizes the translation motion of the joint. However, to

a
use Equations (2.2.3 and 2.2.4), the vector r0 should be represented in the

(X _1, Y“, Z", 1) coordinate system shown in Figure 2.2. 1.

33

2.2.2. Muscle Modeling and Architecture

As proposed in Gordon et al., (1986), the mechanics of muscles are divided into two
major components, 1. the activation and 2. the contraction mechanics. Activation
dynamics is proposed by the work Zajac (1989) to be pararnetrically independent of the
system states (velocities and positions). A detailed discussion of the EMG-to—activation
process and its equivalent signal processing procedure used here was given in Section 2.1.
Utilizing the output of the activation dynamics as the input to the muscle contraction pro-
cess, Zajac (1989) presented a comprehensive development of a non-dimensional muscle
model based on an earlier work done by Gordon et al. (1986). The output of the model
was the muscle force detected at the tendon.

The lumped muscle model proposed by Gordon et al. (1986) shown in Figure 2.2.2,
satisfies four criteria to be suitable for computer simulative studies. The first criteria is the
low order of the mathematical model. Secondly, the model is developed based on muscu-
lar architecture. The third is the model was non-dimensionalized with respect to muscle
and tendon mechanical (static) properties. Last, and most important, the muscle model
output, muscle force, should be consistent with the computed nominal muscle force
obtained using optimization techniques. Static properties of muscle fibers, such as mus-
cle optimal length, optimal muscle force, and tendon slack length, have been intensively
studied. Yamaguchi et al. (1990) surveyed these studies and provided a tabulated sum-
mary of their results for a wide number of human muscles. As has been defined in Chap-
ter 1, the lumped muscle architecture can be divided into two major components, active
and passive elements. The passive element includes both the muscle belly passive ele-

ment and the tendon spring-like element. The active element represents the actin/myosin

34

cross-bridges that are the sight of the activation input (Hatze, 1973).
A typical musculotendon architecture is shown in Figure 2.2.2. In the model shown,
PE is the muscle parallel element, SE is the muscle series element, T is the tendon, and

CE is the muscle active element. The or—angle shown is the muscle pennation angle.

I‘— L,——>i

 

AV

 

 

 

Figure 2.2.2. Schematic sketch of the muscle architectural model developed in
Gordon et al., 1986.
The muscle parallel element is defined as the activation independent elastic component
of the muscle belly. This element has nonlinear exponential properties as shown by Fung
(1967) and Haut and Little (1972). An experimentally based constitutive equation of the

parallel element was developed in Fung ( 1967) that takes the form

35

. (c2 - (Lm - (c3 - L0))/(c3 - L,»

p = c1 ~p0 e (2.2.5)

pe
where p0 and L0 are muscle optimal force and length, Lm is muscle fiber length, ppe is
muscle parallel element force, and c1, c2, c3 are model constants found to be 0.05461,

4.0, and 1.0, respectively. It should be clear that the muscle parallel element will pro-
duce no force for muscle lengths less than the slack length of the muscle. The muscle
series element (SE) is an activation dependent elastic element that was first observed by
the work done by A.H. Hill (1938).

The significance of the series element arises in cases of transient force loading such as
heel strike in running or the takeoff stage in vertical jumping (“Winters and Starks, 1987).
Based on a set of quick release experimental tests, Pierrynowski and Morrison (1985)

developed a series element compliance equation that took the form

Lse/LO = 0.5+(A-x+B-x2+C-x3) (2.2.6)
where
A: 0.2118
B =-0.2265
C= 0.0844
x =pse/Po

and Lse and pse are the series element length and force, respectively.

The constitutive equation for the activation dependent element of the muscle, the con-
tractile element (CE), was first proposed by the work of A.H. Hill (1938; 1950) in the

form of a force velocity hyperbolic equation. This equation relates the contractile element

36

contraction (shortening, Lm < L0) velocity to the force generated through the element. To

account for the lengthening (Lm > L0), Audu and Davy (1985), constructed a bihyper-

bolic equation that can be expressed as

 

r'bl'(1-pce/piso)/(b2+pce/piso) / $1
. pce piso
shortening
V“ = i b 1 / / 133 / (2'27)
_3'( _pce piso) ( ' _pce piso) / >1
. pce piso
. lengthening

where pee is the contractile element force, b1, b2, and b 3 are model constants that takes
different values for different muscles. A. H. Hill (1938) noted that the maximum contrac-

tion velocity in the contractile element Vm occurs when the muscle is at rest and from
Equation (2.2.7), the ratio of b 1/ b2 is equal to Vm . The maximum shortening velocity is
found to be dependent on the muscle optimal (slack) length, L0, and is taken to be 10 L0

sec" (Zajac, 1989). In Equation (2.2.6), the constant b3 was introduced to insure

smooth matching at the point where the two curves, shortening and lengthening, meet.

Values of the different constants in Equation (2.2.7) are tabulated in Audu and Davy

(1985), for selected lower extremities muscles. The isometric muscle force, piso,

appearing in the velocity-force relation given above is taken from the isometric tension-
length curve first observed by Gordon et al. (1966). By investigating the relation between
the isometric force in frog striated muscle at maximum activation, Gordon et al. (1966),
found that the maximum isometric muscle occurs at the muscle optimal (slack) length and
vanishes at above 180 percent and below 58 percent of the muscle slack length. Partial

activation of the muscle will then logically produce less force. Hatze (1977) and later

37

Winters and Starks (1985) suggested that the muscle isometric force-length relation at par-
tial activation is only a scaled version of the relation at full activation provided that the
other characteristics are preserved. Therefore, Hatze (1977) fitted experimental isometric
force-tension data of a fully and partial activated muscle fibers and developed a force-

length relationship at any activation level. Thus the relationship takes the following form

fl

( (we).(4.[Lm_,j

 

 

 

 

-— =4 a +0 -e -— t'00) (223)
pa 1 2 L0
A J 2
where
01 = 0.32
02 = 0.7]
03 = -1.112
04 = 3.722
as = 0.656

p0 is the muscle optimal (maximum) isometric force, and a(t) is the muscle activation.

Finally, the tendon elasticity is assumed to exhibit the same exponential form that the

muscle parallel element exhibits.

A full mathematical description of the lumped model shown in Figure 2.2.2 is given in

Appendix A and takes the general form

p“ = f(Lm,. Vm,.p,a(t)) (2.2.9)

38

where p‘ = g—f , Lm, and Vm, are the musculotendon length and velocity, respectively.

To avoid the complexity and the nonlinearity of the lumped muscle model discussed
above, Khang and Zajac (1989 I, 11) developed a linearization procedure through which a

linear model of the form
p’ = -c ~ p(t) + c - p0 - a(t) (2.2.10)

was constructed. The constant (c) appearing in Equation (2.2.10) above was determined
by the rate of tension of the muscle. To determine the rate (c), Equation (2.2.9) was sim-
ulated using the muscle and tendon static properties assuming that musculotendon length
and velocity were constant at maximum activation (a(t)=1). To illustrate the computation,
Equation (2.2.8) was simulated for the rectus femoris and the hamstring muscles. The
simulation starts with the computation of the musculotendon length, Lm. The muscle was
modeled as a straight line connecting the origin and insertion points. The musculotendon

length can be easily computed from Equation (2.2.2) as follows
1
Lmt = lRO/Il (2.2.11)

A graphical display of the simulation result of Equation (2.2.9) appears in Figure
2.2.3 below. Keep in mind that the linear model of Equation (2.2. 10) is a first order math-
ematical model with a time constant of I/c. Thus, the time constant of the force develop-
ment in the rectus femoris muscle, for example, is the time it takes the muscle force to
reach 63.21 percent of the steady state value, see Figure 2.2.3.

This new linearization procedure used in the present work was adopted for two rea-

sons. The first is based on the nature of the model developed, where every variable

39

(Lm V," p, a(t)) explicitly defined in the model is an implicit function of either the

system states (angles and velocities through the computation of the musculotendon length,
see Equation 2.2.11) and/or the other explicit variables appearing in the model. These
interdependencies make the J acobian based linearization procedure hard to perform in this
case. The second reason depends on the way the force model is used in the context of the
present study. The form given in Equation (2.2. 10) will not be used in its present form but
a sequential transform of the equation instead is used. Also the assumption is made that
the muscle force rate of development from one value to another is constant over the sam-

pling time period (1/1000 of a second.)

40

 

 

 

 

 

 

 

 

 

I r r r r r f 1
o.e - ‘_ fl _, .. - -.
/ ’ a
0.. '- I, —r
"i ’1
a ”
— I
3 0'7 — c-50937 ” -
r
a z
_. .................... .1. ..................................... ..
9" o.e— 5 xi -
\ I ’I ! C - 3e4297
& I I l
v 1 x a
2 , : -
r
r“ i,’ i
l I
E 04r- ’! ! -
' I
E , : : —— Hamstrings
I . .
I I
° / i i
Z o.e- , i i ---- RectnsFemor-is -
I e .
1 I I
,I i i
I I
0.2 ’1' E E
' ! 1
f l l
I I
0.1 "I E E
r l l
I I
0 1 I 1 I1 1 1 1 1
0 0.1 0.2 0.0 0 l 0.‘ 0.. 0.7 0..
Time (seconds

Figure 2.2.3. Simulation display of the muscle nonlinear contraction dynamics
of the rectus femoris and hamstring muscles. The rate of tension
defined in equation (2.2.10) is computed as shown for both muscles.

41

2.2.3. Musculoskeletal Mathematical Model
The general governing equation of motion of the musculoskeletal model can be
expressed through the direct application of Lagrange-Euler formulation. Consider the

general form of Lagrange-Euler equation as

3T av
.687? ) ‘a—qu' 34...; -1: (2.2.12)

where qi represents the i‘h generalized coordinate. In systems that consist of articu-

lated rigid bodies, the human structure as an example, qi may represent either the seg-

mental or joint translational and/or rotational degrees of freedom. In the context of the
present study, it is appropriate to choose those degrees of freedom when modeling the
open loop system that were the subject of measurements during the process of the closed
loop system identification. Appearing on the right hand side of Equation (2.2.12) is the
segmental or the joint applied moments, see Equation (2.2.4) for details.

Body segmental inertial parameters are essential in the formulation of the kinetic (T)
and potential (V) energies expressed in Equation (2.2.12). Segmental masses and the loca-
tions of the center of masses, are usually estimated based on either cadavers studies or
using geometrical and material approximations (Hanavan, 1964; Clauser et al., 1969;
Seirge and Arvikar, 1989).

The derived differential equation of motion can be represented in a matrix form as fol-

lows
1(9) - ('1' = H(q, q) + C(q) + T(q. p) (2.2.13)

where q are the generalized coordinates (outputs) q e R', H (q, q') are the nonlinear

42

Coriolis and Centripetal force vector, H (q, q') e R' , C( q) is the gravity loading force
vector, C(q) e R' , T(q, p) is the generalized force and/or torque vector,
T(q, p) e R' , I (q) is the system inertia matrix, and p is the muscle forces vector

(input), p e R”. The acceleration vector 4 of Equation (2.2.13), can be expressed as

('1' = f(q.4.p) (2.2.14)
where

f (q. q. p) = J"(q) - (H14. 4) + C(q) + T(q. 12)) (2.2.15)

provided that J-l(q) exist. The output feedback estimation procedure proposed here is

based on a family of linearized perturbation equations of Equation (2.2. 14). The family of

linearized equations are parameterized by a set of nominal (operating) points that define

the nominal generalized coordinates 4m , velocities q'in and forces pin. Using Taylor

series expansion of Equation (2.2.14) about the nominal points and neglecting the higher
order terms of the expansion, the associated linearized perturbation acceleration of the
model is

.q' = qu lno8q+qu In-fiq-I-fo In-5] (2.2.16)

where Vt, f I", V q f I , and Vp f In are the Jacobian matrices evaluated at the nonri-
n

nal generalized velocities, coordinates, and muscle forces, respectively, 5q' = q - q" ,

Sq = q — qn , and 8p = p — p" . The smoothness of f (q, q', p) and the assumption

that fiqji«8qj, qut<<5qr and 5pji<<5pj for V IZO and i = 2, 3,..., provides the

43

necessary conditions to neglect the higher order terms in the Taylor expansion. Thus, the
linearization procedure can be used in human activities, such as walking, squatting, etc.,

that do not violate the necessary conditions for the approximation Equation (2.2.16).
The nominal generalized coordinates vector, q" , is obtained experimentally by using

either electrogoniometric systems (Finley and Karpovich, 1964) or photogrammetric

means (Abdel-Aziz and Karara, 1971). Nominal velocities, q', , are computed from
measured generalized coordinates trajectories.
The nominal muscle force vector, pn . Can be computed using the linearized muscle

contraction dynamics Equation (2.2. 10). The activation input, a(t), is computed from the
electromyosignal, measured experimentally, after using the equivalent to activation
dynamics procedure explained in Section 2.1. However, initial muscle forces are essential
to the computation of the nominal forces when using Equation (2.2. 10). Thus, a nonlinear
quasi-static optimization procedure is used to solve for the initial, or time zero, values of
muscle forces. TWO functional quantities need to be specified in order to perform the opti-
mization, namely the objective function and the constraints equations. For the redundant
musculoskeletal system, Crownishield (1978) hypothesized that the redundancy in the
muscle forces spanning a joint is such that the muscles will distribute forces to minimize
muscle energy. In later work, Pandy er al. (1990) proposed a mathematical model to

define the minimum muscle energy hypothesized by Crownishield in the following form

6 = 2 (pi/12,92 (2.2.17)

i=1

where pi is the i’h muscle force and pic is the i'h maximum isometric muscle force. Static

44

equilibrium of a joint may be chosen to be the set of equality constraints by which the
objective function provided in Equation (2.2.17) is minimized. Hence, the form of the
constraints that define static equilibrium can be easily obtained by setting both the gener-
alized coordinate velocities and accelerations equal to zero in Equation (2.2.13) leading to

the following general form
cur.“ = 0)) + T(q,.(r = 0112.0 = 0)) = 0 (2.2.18)
It should be clear when using Equation (2.2.18) that the initial system configuration is
known from experimental data, and the goal of the Optimization is to seek the values of
pn(t = 0) which minimizes (D.
To implement the position feedback estimation procedure defined in the next section,
Equation (2.2.16) is discretized by using the Laplace-transform and the bilinear transform

(see Lewis, 1987). Then, Equation (2.2.16) can be expressed in terms of the backward

delay operator 9, (see Section 2.1) as follows

110(9) . qiq) = 1m) - p(q) (2219)
where
A°(q) = A; - 9‘2”: - g“ +1, (2.2.20)
30(4) = 39944.3? {1+3}; (2.2.21)
and
B? = 610*)“ -B,.* i = 0,1,2 (2.2.22)
A? = (110*)‘1 -A,.* i = 1,2 (2.2.23)

In Equations (2.2.22 and 2.2.23) the matrices are related to the system Jacobians as

45

follows
110* = V, f |n+(2/r_,).vq f ln+(2/ts)2-I, (2.2.24)
Al“ = 2~Vq f |n—2-(2/r,)2.1, (2.2.25)
A2* = V, f ln—(z/r,)-Vq f ln+(2/ts)2-I,. (2.2.26)
30* = 32* = B,*/2 = V, f [n (2.2.27)

where t, is the sampling interval, I, is an (rx r) identity matrix, and (r) is the number of

outputs.

The order of the open loop difference equation given in Equation (2.2.19) is clearly of
a second order. Since the constant output feedback does not change the order of the open
loop difference equations (Fallside, 1974; Lewis, 1987), the order of the closed loop
system (n) to be identified is thus chosen to be 2. Both the closed and open loop systems
should not only have the same order, n, but also should be modeled using the same I/O
structure, ARX model for example. These two conditions are fundamental in the concep-

tual implementation of the position feedback estimation procedure.

2.3 Estimation Of the Position Feedback Structure

An alternative to using the state feedback and the optimal control methods in the study
of the feedback control problem in humans involves the estimation of the feedback struc-
ture. This method is based on the identification of the closed loop system from measured
nominal trajectories of both the input and output data and on the linearized perturbation
equations of a proposed open loop structure.

The block diagram shown in Figure 2.3.1 illustrates the general concept of the feed-

46

back structure estimation proposed here. The matrix [K], shown in Figure 2.3.1, defines
the feedback structure gain matrix that relates the system outputs with the system inputs.
Since we measured only position data in the identification process presented in Section
2.1, the feedback structure will represent only the position feedback control of the open
loop system. However, if both the position and velocity data are available for measure-

ments, the gain matrix [K] can represent the general state (positions and velocities) feed-

 

 

 

 

 

 

back matrix.
I" ———————————————————— 1
Activation ' Generalized
Trajectories : Coordinates
—P [G] _(¥>_> Linearized Perturbed Loo
U Equations in (2%) p > q

 

 

 

 

 

 

 

Feedback Structure
[K1

 

 

 

 

Identified Structure in Equation (2.2.5)

Figure 2.3.1. Schematic illustration of the output feedback
estimation procedure developed in this study
The first step in the development of the output feedback estimation procedure is to
construct a state space realization of both the open loop, Equation (2.2.18), and the
closed loop, Equation (2.2.5), systems. Regardless of the form of the realization chosen,
both systems should be represented using the same state space structure. Utilizing the

block observer form, see Appendix (B), Equation (2.2.5) is written as

47

X(k+ 1) = Ac-X(k)+Bc- U(k) (2.3.1)

q(k) = Cc - X(k) + D". U(k) (2.3.2)
Similarly, the linearized perturbation equations of the open loop system of Equation

(2.2.18) can be written in the block observer form as
X(k + 1) = A” - X(k) +13" . P(k) (2.3.3)
q(k) = C” . X(k) + D" - P(k) (2.3.4)
where X (k) is the state vector, X (k) e R" , q(k) is the output generalized coordi-
nates, q(k) e Rr , P(k) and U (k) are the force and activation input vectors, respec-

tively, and P(k) e Rm ; U (k) e R” . The matrices given in the open loop state space
model are dependent on both the operating points, through the extended linearization,
and the model static and dynamic parameters.
The output feedback law can then be written as follows:
P(k) = G- U(k)—K - q(k) (2.3.5)

Substituting Equation (2.3.5) into the output Equation (2.3.4), leads to
—1
P(k) = [1... + K . 0"] ~{G- U(k) —K - C" - X(k)} (2.3.6)
where I m is an (m x m) identity matrix. Equation (2.3.6) applies only if the inverse of

U... + K . 0"] matrix exists, in other words the determinant of [1m + K - Do] it o . This

condition is defined in Chen (1970) as the well-posedness of the output feedback problem.
As will be seen later and for all practical reasons, this condition is satisfied in the present

study.

48

Substituting both Equations (2.3.5 and 2.3.6) in the state Equations (2.3.3) yields

X(k+l) = {no—13".[1,,,+K-D"]'l ~K-C°}-X(k)+

{3° . [1m + K - 1901’l . G}- U(k) (2.3.7)

By direct comparison of Equations (2.3.7) and Equations (2.3.1), a closed form solu-
tion of the [K] and [G] matrices can be determined as follows (Munro, 1974; Fallside,

1974)
o o "1 o c T o T '1
B-[Im-i-K-D] -K=(A —A)-C° ~(C-C0) (2.3.8)
and
-1
3" - U... + K - D”) - G = BC (2.3.9)
7'

A detailed discussion on the existence of the inverse of C0 . Co appears in Appendix
(C). The necessary and sufficient conditions of the solution of Equation (2.3.8) for [K] are

found in Munro (1974), and Fallside (1974). The solution is also extensively discussed in

Vardulakis (1973). A unique explicit solution of Equation (2.3.8) for the matrix [K]

r —1
depends on the existence of the (8" 3") matrix, see the proof given at the end of this
section. Since the B" matrix is an operating point dependent matrix, or, in other words,

7 -1
an activity dependent matrix, the (Bo BO) may or may not exist for a particular activity

(see later sections). Thus, an approximate least squares based method is considered in the

present study. Hence, a minimum norm 2 solution of Equation (2.3.8) is given as

49

K = 1m.B- 11,—00-81'1 (2.3.10)

where
B = B°*-{(A°—Ac) - COT-(Co- C054} (2.3.11)

and 30* is the generalized inverse of Ba. It is well known in linear algebra (Ortega,

1987) that if Bo is a real (nr x m) matrix, then there exist orthogonal matrices

M = [m1,m2, ...,mnrle R"""”, N = [n1,n2, ...,nmle Rmx'" such that
~ 3" = M-A-NT (2.3.12)
A = [[diag(0'1,0’2,on,1,0,0...)] o] e R”""";Il Snr (2.3.13)

where 11 is the rank of Ba and 31>02> >Onrt >0 are the singular values of the

matrix. Thus the generalized inverse of B" can be represented in terms of M, A, and N

asfollows
30* = N-A'LMT (2.3.14)

A, _ [[diagU/ol, 1/0'2, l/O

tr.» 0’ Own] 6 12"""";1l Snr (2.3.15)
0

It is important to note that the desired closed loop matrix Ac can only be fully attained

if the row rank of B" is full (i.e raw rank(B°) = II , where I1 = nr). For simplicity,

consider the case where D0 = 0. Also let A be defined as follows

50

A =A0-Bo-K- C" (2.3.16)
where A represent the closed loop matrix to be attained, and [K] is the minimum norm 2
based output feedback solution given in Equation (2.3.10) provided that Do is a zero
matrix. Thus, the term [1,, + K001"l in Equations (2.3.9 and 2.3.10) is only an (m x

m) identity matrix. Hence, using the result in Equation (2.3.10) for K and the definition

in Equations (2.3.12 and 2.3.14) for Ba and 80* , respectively, in Equation (2.3.16)

yields
A =Ao—M-A-NT-{NoA*-MT-(A°-Ac) - COT- (C0-C"T)'l}- C" (2.3.17)

In addition to N and M being orthogonal matrices, for the block observer realization

chosen in this study, see Appendix C, the following is true

 

 

—1
(C"-C"T)_l = [[11,] [0]] - [ED = 1, (2.3.13)
and
"11111,," 101
A°-A° = [A2]”, [0] (2.3.19)
..[Anlrx' [0]. nrxm'
where

A, = .43 ..Af, i 1, 2, ...,n (2.3.20)

51

Note that n in the above equations represents the order of the difference equations and
r is the number of outputs. Using the results in Equations (2.3.18-20), the orthogonality

of V matrix, and the definitions in Equations (2.3.13 and 2.3.15), Equation (2.3.17)

reduces to

[Alerr

o I x
A=A —M- I" ’1] [0] -MT- [A2],“ [0] -C”T-C° (2.3.21)
[0] [OIZXIZ] annr

 

 

1A,],x, 101

1.. -annr

where

COToC" = = [[1,1] - [[1,] [0]] = [Ur] [0:] (2.3.22)
[0] [0] [0 annr

and II is the row rank of B” or the number of the non-zero singular values of 3° with

12 = nr—l1 20. Now in the case where I1 = nr, i.e B" is ofa full raw rank, together

with the orthogonality of the M matrix and the result of Equation (2.3.22), Equation

(2.3.21) reduces to

A =A"- [A2]," [0] (2.3.23)

 

 

.J "I“ X [If
It can be easily shown by using the definition of A0 from Appendix (C) and Equation

(2.3.19), that A = Ac. Therefore, the closed loop matrix Ac can be fully attained only

52

if B" matrix is of a full rank. It is also obvious from Equation (2.3.21) that in the case
where 11 < nr the minimum norm 2 solution will only affect and change the II x 11 sub-
matrix of A0 matrix, i.e incomplete attainability of the closed loop matrix.

It was shown from the previous discussion that the row rank of 8” matrix plays a

major roll on the accuracy of estimating the feedback gains. Physically, the number of
independent rows of Ba depends on the type of activity under study. Consider the gen-

eral form of Ba matrix, in terms of system inertia, Coriolis, and gravitational matrices,

given as follows

'5—
B" = __1 (2.3.24)
32
where
— 2
131 = —(A0*/2)-((V6 f ln-(Z/ts) .Ir).30*+30*) (2.3.25)

5-2 = -(A0*) - [[VO f ln—(Z/ts) . VO f In + (2/ts)2 - If} 30* + 30*] (2.3.26)

where A 0*, and, B 0* are given in Equations (2.2.24 and 2.2.27). The Jacobian matrices

V9 f , and, V 9 f are given in terms of system positional and velocity data in C.13

and C. 14. Since V O f is a velocity dependent matrix (see Appendix C for details), it is

clear that it has more effect on the raw rank of Ba in activities that exhibit higher velocity

53

values, i.e ballistic.
Similar argument can also be made for the solution of Equation (2.3.9). Given that the
output feedback problem is well-posed, the solution of Equation (2.3.9) can be expressed

explicitly as follows
G = [1m + K - D0] - 3"” - 3“ (2.3.27)

Recall from Section 2.2 that B" and A0 are functions of the operating point matrices.
Consequently, the feedback matrix [K], found from Equation (2.3. 10), is also a function
of the system operating points, hence the feedback structure is said to be adaptive. Thus,
it is now possible to trace the feedback structure matrix as a function of the operating
point. Because of the large number of the coefficients (elements), r - m , the synthesis of
these coefficients in relation to the inner-connections of the system should be carefully
examined. Several possibilities of input/output connections, (see Equation (2.3.7)), will

be discussed in the next Chapter.

CHAPTER THREE

AN ILLUSTRATIVE EXAMPLE

In this chapter, the proposed position feedback estimation procedure is applied to a
selected human activity, squatting. A detailed discussion of the modeling aspects of the
open loop system (skeletal system) is presented in Section 3.1. Issues concerning both the
selection of the skeletal parameters, masses and moments of inertia, the muscle parame-
ters, muscle static parameters and origin and insertion data are also addressed in this sec-
tion. The modeling of the closed loop system is given in Section 3.2. In this section, the
experimental setup together with the results of input (muscle activation, EMG) and out-
put (segmental motion) data for the squatting test are shown. Results of the identification
procedure are also given in this section. Lastly, the position feedback estimated data are
shown in Section 3.3. Also given in Section 3.3 are interpretations of the data obtained
from the solution of the position feedback problem. Finally, a discussion of the modeling

concerns in the estimation of the output feedback matrix are presented.

3.1 AnthrOpomorphic Musculoskeletal Model “Open LOOp”

A four segment model representing the squatting activity consists of the upper body,
the thigh, the shank, and the foot. Since squatting is predominantly a sagittal plane
motion, each body segment was assumed to have a single degree of freedom. Hence, the
joints, ankle, knee, and hip, were modeled as single frictionless hinge joints. The ana-

tomical locations of the joint centers were selected using the same definitions given by

54

55

Brand et al. (1982) and Hoy et al. (1990). The upper body joint location (the hip joint in
the present study) was located at the center of the acetabulum. The knee joint was located
at the midpoint between the medial and lateral femoral epicondyles. Finally, the ankle
joint was considered to be at the midpoint between the medial and lateral malleoli. All
anthropometric measurements of muscle origins and insertions given in Brand et al.
(1982) and Hoy et al. (1990) are with respect to coordinate systems located at the joints’
centers. Their data were directly used here without any modifications.

The planar linkage proposed assumes a left to right symmetry and is constrained by a
joint located at the toes. Figure 3.1 provides a schematic stick figure of the proposed

model.

Gluteus Maximus V

  
 

Thigh
H .
amstnngs \ Vastr
1‘
I
Gastrocnemius
Shank
Soleus
, Tibialis Anterior
E901

Figure 3.1.1 The sagittal plane linkage model including a set of muscles and
muscle groups representing the squat activity under study

56

In this study, muscles with the same sagittal plane torque function, and which could
not be lumped together, were considered as separate muscle groups. Since surface elec-
tromyography (SEMG) was used to measure the myosignals, we were limited to includ-
ing only superficial muscles in the model. Hence, the total number of muscles and
muscle groups that are the dominant muscles to control the sagittal squat model as well as
being superficial was seven, see Figure 3.1.1. The muscle groups were, from proximal to
distal, gluteus maximus, hamstrings, rectus femoris, vasti muscles, gastrocnemius,
soleus, and tibialis anterior. Similar models were considered in the study of the postural
and the vertical jump control problems by Gordon (1991), Anderson (1992) and others.
TO deal with muscles that have more than one insertion or origin point such as the Tibialis
Anterior muscle, effective origins or insertion points were introduced to best represent
the musculotendon path. That is, a curved musculotendon geometry was approximated
by a set of straight line elements connecting the effective origins and effective insertions.
Values of effective origins and insertions are also provided in Hoy et al. (1990).

The muscle static parameters were modified from those in the literature by comparing
the net muscle isometric moments about a joint with experimental moment data. The out-
come of such analysis affects the simulation results of Equation (2.2.9) for the evaluation
of the rate of muscle force build up (c) appearing in Equation (2.2. 10). Equation (2.2.10)
was used to evaluate the nominal muscle force data, thus it is very important to choose the
appropriate muscle static parameters to obtain an accurate nominal muscle force.

After the muscles’ attachment coordinates are found, the computation of moment and
moment arm values can be found using the method defined in Section 2.2. For the muscle

groups shown in Figure 3.1.1, the origin insertion data are not presented here, but plots of

57

the isometric moments versus the joint angles will be presented later in this section. By
definition, the isometric muscle moment is equal to the product of the muscle moment
arm and the isometric muscle force. As an example of the moment arm calculation, con-
sider the gastrocnemius muscle of the ankle joint. The insertion point of the gastrocne-

mius muscle is at the foot or in the (X1, Y1) coordinate system, see Figure 3.1.2. Also, its
origin is at the femur or the (X 3, Y 3) coordinate system. As shown in Figure 3.1.2., the
muscle spans two joints, the knee and the ankle joints, hence the (m) appearing in Equa-
tion (2.2.3) is equal to two in this case. The transformation matrix T1 is by definition, the
transformation of the (X2, Y2) with respect to the (X1, Y1) coordinate system, thus repre-
senting the ankle joint degree of freedom. Let the ankle joint angle be 0‘, , then for the 2-
D model shown in Figure 3.1.2, the (X2, Y2) to (X1, Y1) coordinate system transformation

matrix T1 takes the form
T1 = [cosea sinea] (3.1.1)

-sin 0,, cos Ga

Similarly, the transformation matrix of the (X 3, Y 3) to (X2, Y2) coordinate system is
defined as T2 which represents the knee rotational degree of freedom. For the proposed

model this matrix takes the form

T2

[cosOk sinfik] (312)

—sin0k c050,c

The joint angles appearing in Equations (3.1.2 and 3.1.3) are related to the segmental

angles shown in Figure 3.1.2 as follows

58

Ga = (02—01)-1t/2 (3.1.3)
0" = (92-63) (3.1.4)
The vector P1 in Equation (2.2.3) is defined as the origin to origin distance from the

(X1, Y1) to the (X2, Y2) coordinate systems. It is clear from the sketch in Figure 3.1.2. that

P1 = 0 because both origins are modeled to be located at the same point. The distance

-> —>
from the (X2, Y2) and (X3, Y 3) origins, P2 , is equal to the shank length. The re and r,
vectors are the origin and insertion position vectors defined from the origins of the (X 1,

Y1) and the (X 3, Y 3) coordinate systems, respectively. Both vectors were found from the

anthropometric data provided by Brand et al. (1981) or Hoy et al. (1990).

59

m4, I4

   

Hip Joint
X4
\
\ 1 X3
Knee Joint v ’
Y1
Ankle Joint

 

Figure 3.1.2 The sagittal plane linkage model showing the segmental
coordinate systems

60

Hoy et al. (1990), developed a procedure by which the muscle isometric force can be
obtained. This method reduces to a solution of a first order nonlinear algebraic equation.

When the musculotendon length, Ln", muscle static parameters, and muscle architecture

are known, the only unknown is the internal variable defined as the muscle fiber length,

Lm, see Figure 2.2.2. Since the goal of the analysis is to compute the isometric muscle
force, the muscle activation is considered at its maximum, i.e a(t) = I, see Section 2.1.
In comparing the tendon force, p , computed first from a tendon model and then com-

puted from the architectural configuration of the muscle, the muscle fiber length can be
obtained. A linear model of the tendon was used when the muscle is in isometric contrac-

tion and is given as

fi = k ' (it-Ls) (1.15)
where p = p/po is the normalized tendon force, 1:, = Lt/Lo is the normalized tendon

length, I; = Ls/Lo is the normalized tendon slack length and, p0 and Lo are the mus-
cle optimal force and length given in Table 1, respectively. The non-dimensional tendon
stiffness, K , was assumed by Zajac (1989) and Hoy et al. (1990) to take the form

K = 375/2, (3.1.6)

The constant appearing in Equation (3.1.6) is related to the tendon elasticity and its
peak isometric stress. More details of the variables shown in Equation (3.1.5) are in Sec-
tion 2.2 and Appendix A. Using the expressions in A.1 and A3 together with Equation

(3.1.5), the normalized tendon force reduces to

61

 

 

a = 37.5/£.~[£..1—£.-L..- J1 ("2"”)2] (3.1.7)

where (to is the resting muscle pennation angle given in Table 1, Z", = Lm/Lo is the

unknown normalized muscle fiber length, and in" = Lmt/Lo is the normalized muscu-

lotendon length as defined in Equation (2.2.11). The tendon force can be computed by
imposing force equilibrium on the muscle architecture shown in Figure 2.2.2. Thus, the

tendon force can be also expressed as follows

 

 

sin “of

3.1.8
1m ( )

i, = (fiiso(£m) + fip¢(zm)) ’ J1 "(

where p,,0(i.,,,) and 13,41...) are the muscle stretch/isometric force relation given in

Equation (2.2.8) and the muscle parallel element constitutive relation given in Equation
(2.2.5), respectively. Due to the fact that, during squatting, the strain in the muscle and
tendon can be large, the muscle passive (parallel) element was included in this analysis.

Hoy et al. (1990) chose not to include the parallel elements.

The isometric tendon force was computed iteratively for every L," by finding the L",
that satisfies both Equation (3.1.8) and Equation (3.1.7). After the isometric muscle force
is computed, by substituting the value of l," in Equation (3.1.8), the muscle isometric

moment can be evaluated for each positional data point by multiplying the moment arm by
the isometric force. Instead of presenting the data for muscle’s origins and insertions and
isometric forces, plots of muscle isometric moments as a function of joint angles are pre-
sented. These plots were compared against experimental and analytical data provided in

the literature (Inman et al., 1980; Marsh et al., 1981; Lindahl et al 1981; Sale et al.,

62

1982; Nemeth, 1983; Hoy et al 1990) to verify the use of both the muscle static parame-

ters given in Table 1 and the muscle attachment coordinates.

 

 

 

 

 

 

 

 

M 1 Optimal Optimal Pennation 1:12;?
use e
Name Force, p0 Length, L0 Aggle, 010 Length, L:
m e rees

(N) ( ) ( g ) (m)
Tibialis 1500 0.070 5.00 0. 1455
Anterior
Soleus 3599 0.034 20.0 0.2376
Gastrocne- 2372 0.072 12.0 0.4180
mius
Vasti 7020- 0.090 10.0 0.1264
Rectus 1344 0.082 5.00 0.3400
Femoris
Harnstrings 3055 0.220 9.00 0.2100
Gluteus 1798 0.180 3.40 0.0090
Maximus

 

 

 

 

 

 

 

Table 1: Musculotendon static properties of the muscles used in the analysis

Although measured muscle static parameters, muscle maximum isometric force,
muscle fiber optimal length, tendon slack length, and muscle resting pennation angle, are
available in the literature (Brand et al., 1982; Hoy et al., 1990; Yamagutchi et al. 1990;
and others), the parameters of some of the muscles used in the present study were modi-
fied. The adjustments were made such that the sum of the computed isometric muscle
forces at a joint are within 10% of the reported experimental isometric joint’s moments
(Inman et al., 1980; Marsh et al., 1981; Lindahl et al 1981; Sale et al., 1982; Nemeth,
1983; Hoy et al 1990). An exact match between the reported experimental data and the

simulation results obtained in this study is impossible due to discrepancies in both the

63

angle definitions and joint locations between both studies. For example, the ankle plantar
flexion moments given in Figure 3.1.3. show that data points measured by Sale et al.
(1982) are in a agreement with the computed moments in plantar and low dorsi fiexion
ankle angles. However, differences were seen at high dorsi fiexion angles. This may be
attributed to the location of the joint center definition which was not explicitly defined in
the work of Sale et al. (1982). Similar observations had been reported in Yamaguchi
(1989) and Hoy et al. (1990). Differences between experimental and simulated data may
also be attributed to the fact that the present study includes only the superficial muscles
where the electrical activity is experimentally easy to obtain. However, experimental
data of joint moments include all muscles contributing to the isometric moments at the
joint. Figure 3.1.4 illustrates the ankle doris flexion moment reported by Marsh et al.
(1981) compared to the simulated Tibialis Anterior muscle’s isometric moment about the
ankle joint. Figure 3.1.4, the largest differential between experimental and analytical data
occurs when the ankle is plantar flexed. In this study, the simulated data of the Tibialis

Anterior muscle will represent the action Of all the ankle dorsi flexors.

 

 

      
 

 

 

 

 

A
E. Total Muscle Moments
6 . _ Gastrocnemius
*5 ace
0 ... - Soleus
E
g 0 Saleetal. 1983 o o
g mo — o
'5!
O
E o
h ~—-~-
5 too — ,” ‘~~
é -_.___ --~-i’-.....
..... I, \‘
0 {D I! "~.‘ \\
a .0 '- ,’ ‘ ‘ ‘-\_ ‘\
< I” ~ ‘ ~"- \\\
I’I’ ~.‘.~'T":‘--~.
Ifl—‘a’” 1 J 1 ‘\~ 1
° -20 —1° 0 ‘3 2° .0
Ankle Angle: Dorsi (+) (Degrees)

Figure 3.1.3. Maximum isometric ankle plantar fiexion moments, with knee in full

extension.

 

I
I
l

 

Tibialis Anterior
Marsh et al. 1981

 

O

O
0
l

 

 

 

&
I

D
I
I

I
I
I

Ankle Dorsi Flexion Moments (N.m)
fl
0

I
I
I

 

 

 

a 1
—20 —IO .0

Ankle Angle: Dorsi (+) (Degrees)

Figure 3.1.4. Maximum isometric ankle dorsi flexion moments

 

 

65

The contribution of the individual muscles to the summed muscle isometric extensor
moments at the knee is shown in Figure 3.1.5. As shown in the figure the simulated total
extensor moments were very similar to the experimental data reported by Lindahle et al.
(1981). At knee flexion angle of about 60 degrees the rectus femoris isometric moment is
about 20% of the vasti isometric moment, which is consistent with data reported by Hoy
et al. (1990). Experimental data of the knee isometric flexion moments reported in Inman
et al. (1980) together with the summed muscle isometric flexion moments are shown in
Figure 3.1.6.

Experimentally, Inman et al. (1980) measured the resisting moments of the hamstring
muscles by applying an external load on the shank. Thus, the gastrocnemius muscle was
not included or loaded during the experiment, and the experimental data is only compara-
ble to the simulated hamstrings moments as shown in Figure 3.1.6.

Isometric hip extensor moments together with the experimental data reported by N em-
ath et al 1983 are plotted in Figure 3.1.7. It is important to note that in the work of Nem-
eth et al. (1983), the hip joint angles were not defined with respect to anatomical bony
landmarks, hence it is the range not the trend of the hip moments, shown in Figure 3.1.7
that is the subject of comparison between the experiments of Nemeth et al. (1983) and
analytical data obtained in the this study. The rectus femoris muscle was considered as
the hip flexor muscle because the other hip flexors’ electrical activity was hard to obtain
using surface EMG. The simulated rectus femoris isometric hip moment computed here is
compared to the Rectus hip moment modeled in Hoy et al. (1990) as shown in Figure

3.1.8.

 

 

 

 

 

 

Total Extensors Moments
no _ _ _ .. Rectus Femoris '1
S. no — _ - _ Vasti Musch ._
E
,3 ...- O Lindahletal. 1981 —
o
E l2.
O
2
g :00
.2
g ..
‘13 ..
8
g 4.
.0
0

 

 

 

 

 

 

   
  
   

 

 

 

... r r ' ' r

1 co - g— Total flexor Moments "
’é‘

_ _ - Hamstrings

z. 140 -
z _ - ... Gutrocnemius
a ..._ ’,..—;—\\ InmanetaLl980 ‘
E p” “Q
g ...,- x _
fl
¢ __ _
.i .0
a
8 .0 '- —
g on — _______ _

8° - .4

0° 2:,

 

 

 

Figure 3.1.6. Maximum isometric knee flexion moments, with hip in full extension.

67

 

 

 

“c U U I? I r fi
Total Extemor Moments
49° F -
.. _ _ Hamstrings
“° ‘ _ - _ Gluteus Maximus ‘
? Namath et al. 1983
one —- -

 

 

 

1.0

Hip Extension Moments (N.m)

 

 

 

1

5° ‘; '° Hifi‘F‘lexion fugue (De'gt’rees) ‘°° "°

Figure 3.1.7. Maximum isometric hip extension moments, with knee in full extension.

 

 

.0 t i t I t r

 

7° - ‘_ Rectus Femoris -

... Hoy et al. 1990 - Modeled

..—

 

 

 

Hip Flexion Moment (N .m)

 

 

 

 

'; 2° ‘16 Hip Fle‘x'ion Angle'llDegrees)‘ 3° ' ’°

Figure 3.1.8. Maximum isometric flexor moment of the rectus femoris muscle about the
hip joint, with knee in full extension.

68

The dynamic equations of motion for the model shown in Figure 3.1.3 were derived
using the Lagrangian formulation discussed in Section 2.2.3. Following the notation given

in Equation (2.2.13), the matrix form of the equations of motion is
J(G)-é=B(®)-92+C(G)+D-M(9,)-P (3.1.9)

where 9 is the vector of segmental angles (output), 6 e R4, 9, is the vector of joint

angles (inter-segmental angles), 9, e R3 , P is the vector of muscle forces listed in the

same sequence as they appear in Table 1, P e R7 , J(G) is the inertia matrix, 8(9) 92 is

T 0
],Drs

a vector describing both Coriolis and centripetal effects with CT = [91 92 93 94

a (4 x 7) matrix that transforms the joint moments into segmental moments, and M(9,) is

the moment arm matrix computed using Equation (2.2.3). Details of Equation (3.1.9) are
given in Appendix (C). Also appearing in the appendix are explicit representations of the
Jacobian matrices, defined in Equation (2.2.16), for the system equation of motion,
Equation (3.1.9). Finally, the open loop input/output equation is obtained by substituting
Equations (C.7), (C.15), and (C.16) in Equations (2.2.14-16).

3.2 Closed Loop Model

3.2.1 Experimental Setup

In this section, a detailed discussion is presented of the experimental setup to measure
segmental angles (system outputs) and to measure muscular activity (system input) on a
human subject performing a squat activity. After completion of informed consent (IRB 93
- 580), a healthy 23 years old female subject was asked to perform at least two squat

cycles per trial. While performing the squat activity, the arms were crossed in front of the

69

upper body and were not considered as separate linkages. Given the subject’s height and
mass, every body segment mass and center of gravity location was computed based on
anthropometric ratios given in Seireg and Arviker (1989). Moments of inertia were com-
puted based on the data provided by Dempster (1955). Table 2 lists the body-segments
parameters for the subject performing the squat activity. The parameters appearing in

Table 2 are defined in Appendix C at the end of this thesis.

 

 

 

 

 

 

 

 

ngznt ma (kg) Lei (m) Li (m) 1, (Kg m2)
Foot 1.599 0.1 16 0.254 0.01 1
Shank 7.760 0.270 0.477 0.030
Thigh 13.340 0.232 0.410 0.118
Upper Body 49.616 0.392 0.392 7.779

 

 

 

 

Table 2: The body-segmental parameters used for the skeletal model

To measure segmental angles, three retro-reflective markers were placed on the foot,
the shank, the thigh, and the thorax. Palpable landmarks on these segments were chosen
as locations of the three markers and were attached to the skin with hypoallergenic double-
sided tape. Two of the markers were placed such that an anatomical body axis can be
mathematically constructed. The anatomical targets used to form the anatomical axis for
each body segment were from proximal to distal, T1 and T10 on the thorax, medial and
lateral femoral epicondyles at the thigh, proximal and distal shank targets at the anterior
crest of the tibial bone, and finally the medial and lateral rear foot on the foot. The third
target was placed in the same anatomical plane non-collinear with the other two. Figure

3.2.1. illustrates the targeting scheme used on the shank body segment. Four 100 Hz infra-

70

red video cameras were placed such that a calibrated space of 170 cm high, 120 cm
length, and 100 cm width was within the visual field of every camera. The definition of
the laboratory coordinate system and the calibration space are shown in Figure 3.2.1. The
data acquisition tracking of the targets was done using a Bioengineering Technology and
System (BTS) Elite system. The trajectories were filtered using a third order Butterworth

low-pass filter.

71

  

Proximal Shank Target

 

 

Posterior
Shank
Target

 

1.7m

 

 

 

1m

 

Figure 3.2.1. An illustrative sketch of the targeting scheme used in the present study
on the shank and the location of the laboratory coordinate system in the calibration space

72

The muscle activity was measured using bipolar surface electrodes with on-site pre-
arnplification. To insure maximum detectability, the electrodes were located approxi-
mately at the muscle motor point. The motor points for the muscles used in this study
were obtained from data provided in Griffin et al., (1982, Appendix A and B), where the
motor points were defined with respect to segmental dimensions. The 8 mm diameter sil-
ver electrodes were located on the approximate location of the muscle motor points with
an inter-electrode fixed distance of 25 mm center to center. The distance chosen here was
recommended in the work of Fuglevand et al. (1992), where both experimental and ana-
lytical models of bipolar electrode detectability was discussed as a function of both elec-
trode size and inter-distance. A larger distance reduces the bipolar electrode detection
strength, since, as discussed in Zipp (1978), the increase of the inter-electrode spacing
distance reduces the band width of the electrode’s transfer function affecting the rejection
of important underlying electrical frequencies. For the agonist muscles, the rectus femo-
ris and the vasti muscles, the distance between the electrodes was chosen to be greater
than 80 mm to insure minimal cross-talk (see Winter et al., 1993).

After establishing the site of the electrode placement, the skin was shaved if neces-
sary, prepared with alcohol, and dried. Double adhesive-backed tape was used to attach
the electrodes to the skin. A small eight-channel junction box was strapped to the sub-
ject’s back. A fiber optic wire was used to transmit the signal to a receiver with a sam-
pling rate of 1000 Hz. The signal was then filtered using a band-pass Butterworth filter
having a frequency band of 5-200 Hz.

3.2.2 Experimental Results

Using the measured targets’ trajectories, segmental (fixed to segment) coordinate sys-

73

tems were computed for each sample time. As an illustration of the calculation, consider
the shank segment shown in Figure 3.2.1. The segmental anatomical axis was chosen to

be the axis running between the distal and proximal shank targets, Z,. The unit vector

corresponding to the anatomical axis is given as

 

 

.. Proximal Shank-Distal Shank

 

 

 

 

 

k5 = > > (3.2.16)
Proximal Shank — Distal Shank
Let i; be a unit vector in the segmental sagittal plane defined as
_ Posterior Shank — Distal Shank (3.2.17)

 

 

 

Posterior Shank — Distal Shank
Then the other unit vectors corresponding to the shank coordinate system (X 5., Y3, Z3)

are computed from

.. i“ xk
j, = L—, .:' (3.2.13)
'18 X ksl
and
is = fsxk; (3.2.19)

The location of each target used in the analysis is measured with respect to the labora-

tory coordinate system (X b YL. 2,) see Figure 3.2.1 and the three unit vectors define the

coordinate transformation from the laboratory to the shank coordinate systems. It should
be noted that for all the body segments modeled in this study, the segmental coordinate

system was chosen to be such that the Z-axis is superior, the X-axis is anterior and the Y-

axis is taken to form a right handed coordinate system. The segmental angles 0,- were

calculated using the joint coordinate method defined first by the work of Grood and Sun-

74

tag (1983). The construction of the joint coordinate system is based on the fact that the
rotational kinematics of a body in a general three dimensional motion can be represented
in terms of three consecutive rotations (Euler, 1748). As described in Grood and Suntag
(1983), the key to the method is to choose one axis from one body segment and the other
axis from the other segment (the laboratory in this case). The third axis is perpendicular to
the other two axes, commonly known as the floating axis or the line of nodes. The seg-
mental X-axis was chosen to be the first axis and the Y—axis of the laboratory as the second
joint axis. For the shank segment, shown in Figure 3.2.1, the construction of the joint

coordinate system is shown in Figure 3.2.2.

jL

Figure 3.2.2. The joint coordinate system of the shank body segment

Where the k —direction is commonly known as the floating axis and is computed as

k = |—.-—— ? (3.2.20)

75

The rotation about the k -axis defines the angle in the (1901:) plane defined as the shank

angle 02. The angle 02 is then computed from the following

92 = n/2_sin 401.1?!) (3.2.21)

Similarly all the other angles were computed for more than two squat cycles and are
shown in Figure 3.2.3. The angles plotted are the same angles defined in Figure 3.1.2,

hence the thigh angle 03 at full extension takes a value of 90° as shown in the figure.

76

 

‘
o
O

Segmental Angles (Dew)

T

 

'D'unkAngle 04
mgr-Angle 0:
Shank Angle 02
" ‘- ." s
l . .1 \
x. I 1:. I \
x \
\. .\ \
\. 1 .‘
x x. -‘ r-
\ \. \ .l.
x 1' \_ ~- I _
‘ ~‘ ‘ FullExtemion

 

 

 

Time (mseconds)

Figure 3.2.3. Segmental angles

77

The muscle activation curve (system input) was computed using the procedure dis-
cussed in Section 2.1. The measured electromyographic signal was rectified, modulated
and then normalized to represent the muscle activation dynamics. Figures 3.2.4-7 show

the processed experimental muscle activity for the group of muscles included in the

present work, see Table l.

 

 

 

 

 

 

 

 

 

I I T
I . u-
3 " Gluteus Maximus
1 - t 0‘. ‘
r ’ 1
1 ' 1
| I f .
. 1 1 .
. 1 r :
3‘ L . - I
o’ . I o I —
23:! r. ' " l '
l 1
1 | . ‘ 1 I
< 1 1 . ‘ 1 '
' ‘ o . 1 ‘
0.. r- I I I I I —‘
' I I ‘ . \
I 1 I ' I
I . I ‘ : 1 |
1 I 1 ‘ l ‘1
’ | ‘ I ‘ I
0,4 t— ' I '- I \ Ll
" 1 ' 1 1 \ | \ r
= ’ .1 r " ' ' i 1 '
Z ' 1 I 1 l ' o \ t.
' I a ’ ' I ‘
' I s I ‘ ' I ‘ I
.. ' 1 I \ l I l 1 ‘I
on F-\ K ,I ‘ '1 7" ' 1' t
I I I ‘ ‘ I ‘ I
‘-..‘ ‘, \ ,’ ', \ I,
v \_” \' . 7“ ’
J
o l 4 I l l I I
o 1000 none .ooo 40"! coco coco 7000
Time (mseconds)

Figure 3.2.4. Muscle activity plots of hamstrings and gluteus maximus for three squat
cycles

 

78

 

d
-

 

 

Vasti Muscles
Rectus Femoris

 

 

-
I
I
I
I
I
I
i
I
I
\
I
'I
I
I
I
l

1‘
I

 

 

 

1
I
I
r
r
1
1
o
1
a
a
r
l
1
1
r
I
I
1
I
I
a
I
l

0000000

I‘II‘I'
III
’I
----
o
'v
‘II

'II III

P p
._ - . . a ..
1 0 0 0

.353 E1232

 

Time (mseconds)

 

Figure 3.2.5 Muscle activity plots of rectus femoris and the vasti muscles for three squat

cycles

 

 

 

 

mius

G

 

 

.aL-

 

 

_

8000 0000 50M 0000 0000 7000

 

 

:23 $52

up

Time (mseconds)

Figure 3.2.6. Muscle activity plots of the gastrocenimus muscle for three squat cycles

79

 

 

 

 

 

 

 

 

    

 

 

 

I .2 — —'
Soleus
1 ~ ’3‘ _ _ _ Tibialis Anterior -
’ s
=' ‘1
5:- i =
O; 0.0 -! I:
3:3 ' 1
g : I. '4‘ I"
4 5 t 3 x 5 x ‘
0 0 1 : . a ’r _
E l i l 1' 1 o“
o . ' I ‘ '. ' I
I: 1 i l '. r '. 5
E . .. ' z .' -. : z .-
° ‘ ' z : t ;
Z ‘- :’ 1 .’ '- :
I 1 o '. I
l I .1 1' " 1.
0 .I I. 'I ‘I {‘7‘ 3 I :
'1. ’1' ‘. ,- ‘o ‘g ‘1‘ ,’ :
7“ I 7‘ I ‘ I, ‘ ‘\ ’Il ‘\ I
7-“~” a 1 \’ n ...-.7 a ‘\ 1
0° 1000 8000 0000 ‘00” 0&0 00” 7000
Time (mseconds)
Figure 3.2.7. Muscle activity plots of the soleus and the tibialis anterior muscles for three
squat cycles

3.2.3 Autoregressive Closed Loop Estimation
In this section the results of the identification of the closed loop system are given. Fol-
lowing the procedure discussed in Section 2.1, an Autoregressive (ARX) squatting model
Was constructed by utilizing the input/output data shown in Figures 3.2.4-7 and 3.2.3. To
achieve a robust estimation, a static optimization procedure was constructed. The initial
gu'53Sses of the parameters were chosen such that the closed loop state space model (Equa-
ti0‘18 2.3.1-2) is stable, i.e the eigenvalues lie in the unit circle (Strejc, 1981). This
assumption is consistent with the fact that the observed experimental data represent a real
Subject. It should be clear that this static optimization procedure is independent of the

physical system itself but is dependent on the state space model chosen to describe it.

80

Thus, the method proposed was used initially to get a stable starting point for the predic-
tive error estimation procedure.

The estimated closed loop model, Equation (2.3.1), was then simulated using the
measured muscle signals as inputs and the outputs of the simulation were compared to the
measured segmental data as shown in Figure 3.2.8. A collective effect of all the estimated
angles are also presented on the four degrees of freedom stick model shown in Figure

3 -2.9. In both figures the solid lines represent the measured data while the dashed lines
represent the simulated results. A small time shift between the measured and simulated
data during the first squat cycle (compare the dashed and solid lines of the stick figures of
Figure 3.2.9) is noted. This time shift did not occur during the second squat cycle. The
variable time shift is attributed to the fact that the identification used in this study was of
an off-line non-recursive nature. One approach to eliminate this problem is to use an on-
line recursive identification procedure to estimate the system model. However, the avail-
ability of computer codes for on-line recursive identifier for a large MIMO ARX does not
Currently exist. Thus, the results of the identification procedure given in this section are

Considered sufficient in the sense of minimizing a robust estimated error criterion (Ljung,

r988)

Foot Angle 01

T111211 Angle 03 .

81

 

15"

-IO "

-15'

 

 

 

 

O 2000 soon

 

 

 

 

Time (mseconds)

Figure 3.2.8. Measured segmental angles (solid lines) and simulated

‘
'

Shank Angle 02

 

‘ ‘ ‘ ‘
V

- ‘
V

 

 

 

 

 

 

 

 

 

Time (mseconds)

soon

results according to the MIMO ARX model (dashed lines)

82

.33 Baa—25m 8a 3:: vogue 05 .83 none:— 35308 2: 2a 8:: 38 05.
53:8 3:3 5 3.258 been Emu Eocene .«e mop—woe .58 05 .«o manna xoem .m.m.m oSwE

 

83

3.3 Position Feedback Solution

3.3.1 Validity of the Solution of the Position Feedback Problem
Recall that the position feedback matrix [K] was computed based on a minimum norm
2 approximation criterion, see Equations (2.3.10 and 2.3.11). The approximate solution

of the matrix [K] is found from the solution of Equation (2.3.10) as follows, see Section

2-3,

K = 4,3 - [I,—1D"-B]'1 (2.3.10)

To test the validity of the above solution let us first reconstruct Equation (2.3.7) by

substituting back the values of [K] computed above in Equation (2.3.7) as follows

X(k+1) = {iv—13"-[1,,,+K.D"]’l -K.C"}-x(k)+

{30. [1,, + K . 1201'l - G}- U(k) (3.3.1)

where X(k) is defined as the reconstructed state vector computed at a step k when using

the solution of Equation(2.3.10) for the matrix [K]. Theoretically, for Equation (3.3.11)

‘0 represent the closed loop system defined in Equation (2.3.1), the terms
a —1 ..1
A ~B°-[Im+K.D°] -K-C° and {Bo-[Im+K-D°] o} shouldbeequal to

Ac and Bc , respectively. However, the above statement is only true when the Bo matrix
is of a full raw rank, see Section 2.3 for details. The reconstructed state Equation (3.3.1)
and the closed loop state Equation (2.3.1) were simulated using the same set of inputs and
initial conditions and were compared as a validation of the estimation of the position feed-

back matrix [K]. The validation procedure is presented in Figure 3.3.1.

84

 

 

 

 

A The Value of The State
at time 112
25‘
c Q )0
Given Inputs and \ ° .9 0.90:5 T
Nominal States Data 2“ /
3 /
/
Us; / / I ’ I ' T
/ /
5‘ , ’ \Y 3 9 .
.2 // / ’ \ 133°“0 5
a o .,
//
/ I
I : xi
. E I I >
t1 t2

Time

Figure 3.3.1. Schematic representation of the validation procedure

Let ;i(k) represent an arbitrary state computed from the simulation of Equation

(3 .3.1) and let xi(k) represent the corresponding state computed from the simulated

Closed loop system of Equation (2.3.1). The percentage relative error in computing the

reconstructed state ;i(k) is defined as

15.06) - xi(k)

Xi ( k) X 100 (3.3.2)

 

Relative Percentage Error =

85

Both xi(k) and 15,00 are functions of operating points, thus the relative percentage

error is also dependent on the operating point. As an example, the relative percentage
error for two of the states, the trunk and shank angles, are plotted in Figures 3.3.2 and
3.3.3, respectively. From the figures, the percentage error in both the shank and trunk
angle computations range from zero percent to less than 10%. The maximum error of
about 9% in the shank angle simulation was observed to occur at 200 milliseconds, as
shown in Figure 3.3.2. Small as it is, this maximum value occurs only for 50 milliseconds
and the error then drops back to its lower value of about 2-3%. The 9% error corresponds
to about 8 degrees error at a point where the actual shank angle is around 100 degrees, see
Figure 3.2.3. This maximum range of error falls within the standard deviation across trails
in a normal subject’s squat activity (Bemis, 1992). The same conclusions can also be

made for the trunk angle. As it was stated in section 2.3, the source of the error is the fact

the for the simulated squat activity, the row rank of B" was not full.

86

¢
0

 

O I
l l
l
{\l\.
l l

b
1

Relative Percentage Error

   

 

 

 

 

-10 l 1*

‘0 IOOnme ( 1 6° ands) 2W 0 .00

 

Figure 3.3.2. Relative percentage error plot for the computation of the shank angle.

10

 

 

I

b

r
l

I
O
T

Relative Percentage Error
.3
l

 

 

l
O
l

 

 

I

d
O
p

 

A 1 1 A
‘0 100 1 ‘0 zoo 2‘0 .00

Time (mseconds)

Figure 3.3.3. Relative percentage error plot for the computation of the trunk angle.

 

87

3.3.2. Interpretation of Some of the Position Feedback Elements

Since the B appearing on the right hand side of Equation (2.3.10) is a function of the
system Operating points, defined by a set of nominal muscle forces and motion trajecto-
ries, the matrix [K] computed using Equation (2.3.10) is also a function of the operating
points. Thus for the present application, the 7x 4 elements of matrix [K], 4-outputs and
7—inputs, are expected to change as the operating points change. The matrix [K] was

defined in Equation (2.3.11)
P(k) = G- U(k)—K - q(k) (2.3.11)
represents the contribution of the output vector q(k) to the muscle force vector P(k) .

In other words, the element (_kij)’ for the present example i = 1, 2, ..., 7 and

j = l, ..., 4 , of the matrix [-K] represents the feedback gain relating the jth output of the
system, 61-, to the in. input, pi. For example, the first row of the [-K] matrix, (—k1j)

j = 1, ..., 4, represents the feedback gains or components from the segmental angles

91, 02, 93, and 94 to the tibialis anterior muscle force. A positive value of the (—k,-j) ele-

ment implies a positive contribution of the 1M motion on the force produced in the i M
muscle, see Equation (2.3.11). In other words, the joint motion is demanding more mus-

cle force for the proposed open loop system behavior to track that of the identified closed
loop structure. Conversely, a negative value of the (~kij) element indicates that a partic-
ular motion is requiring less force generation. These element can be represented in terms

of feedback gains relating joint angles to muscle forces since joint and segmental angles

are analytically related, see Appendix C.

88

To demonstrate the adaptability of the elements of the feedback matrix [K] to the
changes in operating points, consider the one squat cycle shown in Figures 3.3.4a and
3.3.4b. The figures show the joint angles and the associated stick figure for one squat
cycle. As an illustration only, consider the feedback coefficients mapping the three joint
angles, ankle, knee, and hip, to the tibialis anterior muscle, as shown in Figure 3.3.5. It
is clear that the feedback gains are not constant throughout the range of time considered.
The values of the coefficients change depending on the configuration, the configurational
speed, and the muscle nominal forces in the system. It is well known that structurally the
tibialis anterior muscle is of no relation to either the knee or the hip joints, however Fig-
ure 3.3.5 clearly shows the motion of this joint is related to the force developed in the tib-
ialis anterior muscle. This observation is consistent with the hypothesis that states that all
the muscles involved in a particular activity are related somehow through a global feed-
back structure. This phenomenon is also clear when considering the feedback gain plots
from the knee and hip joints to the soleus muscle and the ankle and hip joints to the vastus
muscle as seen in Figures 3.3.6 and 3.3.7, respectively. The positive or negative contribu-

tion of these motions, positive and negative contribution is simply related to the sign of

the (-k,.j) element, are dependent on the operating point.

 

Joint Angles (Degrees)

Figure 3.3.4. The stick figure and the associated joint motions for the one squat cycle.

89

__
k
L
—

 

 

 

 

J l_.J‘L_' _

 

_. fl

 

 

 

 

 

 

 

 

 

 

1.0 - r ' -*
ExtensionStarts ---- AnkleAngle-(Dorsi+)
...... _._._. KneeAngIe-(Flexion+) -
#K
,I' -.\ ........ Hip Angle - (Flexion +)
..l’ ’ al— p‘. a 3’1
.1”. 1' "\_ 1' I: 0‘
.I ‘ I. '
‘.° b :{’ ‘ ...". -‘
"I \.‘ ...;o'
’ "o ....i
\. 3'1
we — \. .w’ -‘
s. .3;
‘_ o s."
\ .0 0.
o ‘0 ... 3
. -- \. i -l
, J
\. I. f
\ °,
. . §
a -. End of Externion ..I
.6 '— \. '. ,0 -‘
\. ‘ o'
‘. .3
‘. .0. a
- — - — _ - - ~ - \ .‘ 3 '
40 = ~ ~ ‘. - :- .. - -
\ $ 0: ’ ’
x .‘ 3‘ 'i' I
\ °‘ ‘. 0" I '
\ . $ '0 /
‘ ‘ “ l’
.0 — \ “ “ I I I —I
‘ c .
s \ ’ I
‘ s 's .’o ’
~ .u. _4 ...
o l L l l A
o ‘00 1 O” _ 1 C” .0” “00

Tim (mseconds)
(a)

 

"— —1

Extension Starts End of Extension

10— -

 

 

 

 

—10'

 

Ankle Angle Feedback "
- ._ - Knee Angle Feedback

- - - 'IhmkAngleFeedback _

1 1
‘0 1 00 1 .0 ace “0 000

Time (mseconds)

 

 

 

Feedback Coefficient to the Tibialis Anterior Muscle

 

 

Figure 3.3.5 Feedback gain plot from the ankle, knee, and hip motions to the
tibialis anterior muscle.

 

‘0 - Extension Starts End of Extension ~

 

 

 

 

 

 

 

 

Feedback Coefficient to the Soleus Muscle

 

 

 

-.. L- Knee Angle Feedback -
-.. _ _ _ 'IhmkAngle Feedback _
.3 "E. . i. 23.. ac‘c .ee
Time (mseconds)

Figure 3.3.6 Feedback gain plot from the knee and the hip motions to the soleus muscle

91

 

8
r
1

Extension Starts End of Extension

‘
8
l

l

8
l
L

O
l
l

l
8
I
I
\

 

 

 

I

i
8
l

 

__ Ankle Angle Feedback
.. _ - Think Angle Feedback ‘

I
e.

 

 

 

 

Feedback Coefficient to the Vastus Muscle

 

a a a
.0 100 1 .0 zoo 8‘0 000

Time (mseconds)
Figure 3.3.7 Feedback gain plot from the ankle and the hip motions to the
vastus muscle.

Although the hypothesis of internal coordination has been proposed by different
researchers in the fields of neurosciences and human dynamics, the characterization of
such interaction has not been fully investigated. This is largely due to the fact that a gen-
eral human activity exhibits nonlinear dynarnics that are usually analytically complex.
The use of extended linearization to address the nonlinearity of the system dynamics with
algebraic control theory in the present work can be considered as a first step in the direc-
tion of the quantization of the feedback interaction (Dhaher, 1996).

Let us next examine the effect of the individual joint motion on the feedback gain of
the muscle(s) spanning a particular joint. To demonstrate this concept, consider first the

ankle joint. The graph shown in Figure 3.3.8 is a cross-plot of the ankle angle versus the

92

feedback gain associated with the soleus muscle for the squat cycle shown in Figure 3.3.4.
As the ankle angle increases in the dorsal direction, the soleus muscle feedback gain

increases towards the positive value but remains negative all the way up to an ankle angle

of about 30°. Since the soleus muscle is considered a plantar flexor at the ankle joint, the
trend shown in Figure 3.3.8 indicates that the feedback structure is trying to give the

soleus muscle the signal to facilitate the dorsi flexion process. However, as the dorsi flex-

ion angle increases beyond the 30° mark of ankle dorsi flexion, the feedback gain
switches to a positive value to demand an increase of the muscle force to provide both the
braking effect on the excessive dorsi flexion and maintain the dorsi angle ankle of the open
loop model consistent with the measured ankle angle of the closed loop model. On the
other hand, the Tibialis Anterior muscle is receiving more and more negative feedback
gains in order to facilitate the counter action of the soleus muscle as discussed above, see

Figure 3.3.9.

93

 

10"

 

Feedback Gain to the Soleus Muscle
O

 

 

 

 

-zo - -
-u - ..
-‘O 1 1 1 1 1 1
o to 20 so 40 so so
Ankle Angle: Dorsi(+) (Degrees)

Figure 3.3.8 Ankle angle feedback gain plot to the soleus muscle for a squat activity

94

 

'0 I I I I I I

Ext_and

 

.4-

Feedback Gain to the Tibialis Anterior Muscle

-‘r—

 

 

 

-10 I J 1 I 1 1
O 10 20 M 40 50 60

Ankle Angle: Dorsi(+) (Degrees)

 

Figure 3.3.9 Ankle angle feedback gain plot to the tibialis anterior muscle for a squat
activity

211'

95

Similar trends were also observed at the knee joint. As the knee angle increased, the
feedback gain to the Rectus muscle increases to maintain its negative value to facilitate
more and more knee fiexion as shown in Figure 3.3.10. However, as seen in Figure

3.3.10, this trend does not appear to continue throughout the range of knee flexion. At

about a knee angle of approximately 50° the rectus femoris starts receiving positive feed-
back gains to increase the muscle force to counteract the gravity and dynamic effects and

avoid a total collapse at the knee joint. This is clear from the rapid increase of the feed—

back gain values at large knee angles around 120°. Unlike the Rectus muscle, the ham-
strings are acting in the same manner as the Tibialis Anterior muscle at the ankle joint.
The only major difference is that the hamstrings muscle is required to reduce the amount
of force buildup by getting large values of negative feedback gains at large knee fiexion

angles in order to reduce the counter action against the rectus femoris muscle, see Figure

3.3.11.

96

 

 

 

 

 

 

 

 

l I I I l I I F l
70- .-
eo- - !
2
a ..- -
2
g ”I— —1 h
5
[it
5. so- -
8
a:
g 20!- -
3
'5'
.8
g 0
In
-10- -
-mr- -
_” l l l l l I l l l
-20 O 20 ‘0 GO .0 100 120 1‘0 160 1.0

Knee Angle: Flexion(+) (Degrees)

Figure 3.3.10 Knee angle feedback gain plot to the rectus femoris muscle for a squat
activity

97

 

 

 

 

 

 

o I I I I I I I I I
-50 - -
-Ioo - ‘5._ ~
0
g g 2'- Laat Poht
' O
450 - ° -
“I O?)
00
0
° 0
o o
5 -200 - 0° ‘
3 °8
5 °3
e °°
x 450 - :3 -
o
3 00
g co
k, 400 - E -
O
O
% ExLahrt
_‘so 1 1 I l L l l L 1
-20 O 20 40 60 It! 100 120 H0 160 190

Knee Angle: Flexion (+) (Degrees)

Figure 3.3.11 Knee angle feedback gain plot to the hamstrings muscle for a squat
activity

98

The last joint to examine is the hip/trunk joint. As shown in Figure 3.3.12, the feed-

back gain to the gluteus maximus muscle is negative for all hip angles less than 100°.
However, the feedback gain switches in sign after the angle that forces the gluteus maxi-
mus to produce more muscle force to control the large hip/trunk flexion angle. Although,
both the hamstrings and the gluteus maximus muscles are considered to be hip extensors,
the feedback gain plots does not show a similar positive/negative synergy. While the feed-
back gain from the hip joint motion to the gluteus maximus muscle switches in sign, the
hamstring feedback gain remains positive throughout the squat cycle considered, see Fig-
ure 3.3.13. This may be attributed to the fact the hamstring muscle is considered the
major hip extensor muscle and is considered to be active in a positive sense for the whole

rang of hip flexion angles. Also, note the rapid increase of the positive feedback structure

to the hamstrings muscle as the hip joint angle passes the 140° mark. Finally, it is very
important to not that although the hamstring muscle is getting negative feedback gain from
the knee angle, it is receiving positive feedback from the hip angle, see Figures 3.3.11 and
3.3.12. This is attributed to the fact that the hamstrings acts as an action generator at the

hip joint and a facilitator at the knee joint in the extension mode of the squat activity.

 

159I I I I I I I I I I I

IN?"

 

 

Feedback Gain to the Gluteus Maximus Muscle
8
I

 

 

 

 

ExLend
_50 1 l I l I l l L l l I
-20 o 20 40 so an roe 120 no we no

'h'unk‘Angle: Flexion '(+) (Degrees)

Figure 3.3.12 Hip angle feedback gain plot to the gluteus maximus muscle for a squat
activity

100

 

3mI I I I I I I I I I I

250- -4

- to
8 8
I I

‘
8
I

Feedback gain to the Hamstrings Muscle

Ext_and

 

 

 

1 1 1 l 1 l l 1 I 1 J
-20 0 2O ‘0 60 U I M 120 140 I“ 1.0

Trunk Angle: Flexion (+) (Degrees)

 

Figure 3.3.13 Hip angle feedback gain plot to the hamstrings muscle for a squat activity

attiv
mu

bi-al

join
mus

that

gal

all:

M

h:

101

To further analyze the position feedback synergy in the musculoskeletal in squat
activity, let us consider the position feedback cross plots of the agonist and antagonist
muscles included in the model. Structurally, muscles are classified as a uni-articular or a
bi-articular, where the latter means that the muscle spans two joints at the same time. The
rectus femoris muscle is an excellent example because it spans both the hip and the knee
joints. A uni-articular muscle is a muscle that structurally spans a single joint, the soleus
muscle. When comparing the feedback gains of two antagonistic uni-articular muscles
that are associated with the joint angle both muscles span, a linear coupling was observed.
This is clear in the cross plot of the gain associated with the ankle angle to the soleus and
Tibialis Anterior muscles shown in Figure 3.3.14. Since negative and positive feedback
gains are representations of demands by the system for less and more muscle force gener-
ation, respectively, the positive/negative synergy across the uni-articular antagonistic
muscles in squat activity appears to be justifiable. If one of the antagonistic muscles is
creating an action, braking or driving the opposite muscle should facilitate that action.
This is generally true for the uni-articular antagonistic muscles, however, it appears that
in the case of a bi-articular antagonistic muscles the positive/negative synergy does not
necessary hold. As an example, consider the cross plot given in Figure 3.3.15, where the
hamstring feedback gains associated with the knee joint are plotted against the rectus fem-
oris position feedback gains. It is obvious that the bi-articular nature of both muscles not
only produces a nonlinear form of coupling, but also loses the positive/negative synergy
observed in the uni-articular case. Although the feedback gain to the rectus femoris mus-
cle has switched in sign during the squat cycle, the corresponding hamstring feedback

gain maintained the same sign. This trend may be attributed to the fact that unlike the uni-

 

102

articular antagonistic muscles, the bi-articular muscles are supposed to accommodate the

motion of the other joint hence the feedback structure is more complex.

 

 

 

 

 

 

I I I I I I I fi I
Q- ..
3
g zo- _
2
3 o
0 to- ..
5
3
3 a
en
E
o
a o
< “x
o
e O
E 40* 0° -
6::
CD
'5 -2o- -
a
i
In
-a- -
l I l l l l A l 4
-I -O -4 -2 O 2 4 6 0 10

Feedback Gain from the Ankle Angle to the Tibialis Anterior Muscle

Figure 3.3.14 Position feedback gain associated with the ankle angle cross plot of the
soleus and tibialis anterior muscles for a squat activity

 

 

103

 

 

 

 

 

 

2
I 4°“ '
81
ea -1” "' -1
fl
E.
a: 45° _ o oo _
0 o O
'5 o 00
3 0 o
'3 -200 - D on "
I: a 00
<1 0 on
8 O o
g _2“ F- D on -I
O
0 o
o
«.5 o o
E o
o

5 g o
3 -360 - 0 ° -
O O O
2 ° 0°
3 £3 Ext_atart

-400 ~ -
In

460 - -

-‘” l 1 l J 1 J J l l

-30 -20 -IO 0 IO 20 U 40 50 80 70

Feedback Gain from the Knee Angl'e to the Rectus Femoris Muscle

Figure 3.3.15 Position feedback gain associated with the knee angle cross plot of the
hamstring and rectus femoris muscles for a squat activity

 

104

The next set of muscles is the group of muscles that are agonistic, i.e muscles that
have similar actions at a joint, for example, the hamstrings and the gluteus maximus mus-
cle at the hip joint. Since a uni-articular group of agonist muscles does not exist in the
model proposed here, we will only consider the case of agonist muscles regardless of
their structural relation to the joints. The nonlinear coupling observed in the case of the
bi-articular antagonistic muscle is also present in the case of the feedback gain cross plots
of the agonist muscles when at least one of the muscles is a bi-articular muscle, see Figure
3.3.16. This nonlinearity is attributed to the fact that one of the agonist muscles is bi-artic-
ular, the hamstrings. While the uni-articular muscle feedback gain changes sign, see
Figure 3.3.16, the hamstring antagonistic muscle maintained the same feedback gain sign
throughout the squat cycle. This trend was also observed in Figure 3.3.15. However,
unlike the antagonistic case, the agonist muscles act together to generate an action,
ecentric or concentric. The cross plot given in Figure 3.3.16 supports that definition where
it can be observed that the position feedback synergy of the gluteus maximus and the ham-
string muscle appears to exhibit a positive/positive feedback synergy through the exten-

sion phase of the squat cycle.

105

 

 

 

 

 

 

.9. I I I I I
§ 14o - -
2
3 & Ext_aurt
E 120 - ° -
.5 0
2 .59
B O
3 IN '" so .I
.2
U
o
5 .. - 8° _
3 Q
0 Q
.5.“ 6
1: 03
< co - OB -
a o.
E o
0 ‘0 " 0‘? 1
2 or
E
:5 w -
CD
.I
o
_g o
E

-20 - -

Ext_and
_‘o 1 L L 1 1
O 50 100 ISO 200 250 300

Feedback Gain from the 'fi'unk Angle to the Hamstring Muscle

Figure 3.3.16 Position feedback gain associated with the hip angle cross plot of the glu-
teus maximus and the hamstring muscles for a squat activity

 

CHAPTER FOUR

CONCLUSIONS

This study described an alternate yet physically justifiable approach to computing the
feedback gains. This method is based on the idea that any human activity can be
characterized by an input/output mathematical model. The model inputs and outputs
signify measured muscular activity and nominal motion trajectories. This model
becomes the structure that includes both human dynamics and control, and is a physical
representation of closed loop system dynamics. To obtain such a mathematical model, a
detailed discussion of the parametric identification method was developed. The
development was restricted to a general input/output human dynamic model.

The development of an open loop system is a major component in the study
concerning the feedback problem. An anthrOpomorphic musculoskeletal model (open
loop), where muscles are considered as moment generators, was developed. Since the
open loop model that was developed is nonlinear, an extended linearization method was
used to overcome nonlinearity. The Jacobian-based linearization procedure reduces the
nonlinear model to a family of linearized systems, parameterized by the same set of

measured nominal operating points. Also the generalized forms relating Jacobian
matrices with open loop system inertia, Coriolis, and gravitational matrices are derived.

A key step in the procedure is to represent both the estimated closed loop system and

106

107

the proposed open loop models in the same state space structure representation. The
block observer form was used as a standard state space realization of closed and open
loop systems. A detailed development of an explicit solution of the feedback gain matrix
was presented. Also the conditions for good feedback estimations were addressed. The
utility of the proposed method was demonstrated by applying it to the squat activity.
Quantitative results of the positional feedback gains were presented.

The accuracy in estimating the position feedback matrix is dependent on the type of
activity under study. For more ballistic movements a more accurate estimation is
expected. Since the squat activity is slow, its estimation of position feedback was not
perfect in the sense that all the closed loop model matrices were not fully attained by the
proposed output feedback law. To test how accurate the estimation was, the open loop
system combined with estimated position feedback was first simulated. Then, the result
of the simulation was compared with nominal motion trajectories obtained from the
estimated closed loop model. Although the squat activity studied was slow, the error in
the feedback gain estimation was low.

Compensatory functions of muscles were characterized by observing the adaptability
of the feedback gain values relating a specific joint motion to the force generation in
different muscles in the structure. Results showed that activity of a muscle is affected by
feedback from angular motions of all the joints in the structure. When comparing
feedback coefficients of two antagonistic uni-articular muscles, the coupling was

aPPI'OXimatcly linear. Positive feedback to the driving muscle was accompanied by a
negative feedback to the antagonistic muscles. For bi-articular antagonistic muscles, the

feedback coupling was more complex. The results also showed that a positive feedback

 

108

to the driving muscle at a joint may or may not be accompanied by a negative feedback
of the antagonistic muscle. Positive/positive or positive/negative feedback synergy in
antagonistic bi-articular muscles was shown to depend on many parameters that include
the motions of one or more joints.

Using the proposed method, the feedback structure can be estimated for a general
nonlinear dynamic structure. It can compute the feedback gains without resorting to any
assumptions over the dynamic behavior of the closed loop system. In addition, the
method provided quantifiable data that is valuable information concerning feedback
synergy in human structure. Improvements on this method may include incorporating
recursive (on-line) methods in estimating the closed loop model.

An area of future research may be to investigate the potential use of this method in
solving the complex orthopaedic problems of predictive surgeries. Such problems
involve the prediction of the motion of a human structure after undergoing a muscle
transfer surgery. It is well known (Vukobratrovic er al., 1990) that the open loop
human structure is inherently unstable in the absence of any feedback control. It has
been shown in this study that the estimated feedback gain, together with the open loop
system, results in a collective system that simulates the real activity, the closed loop
stable system. Therefore, surgery can now be simulated on a stable system that consists
of the open loop model and the estimated feedback gain structure. Questions pertaining
to the effect of surgeries on muscle activation and position feedback gains should be
addressed. Studies appearing in the literature showed that muscle activation does not
vary with the transfer surgery (Waters et al., 1982, as an example). Another issue of

concern is the sensitivity of the pre-surgery estimated feedback structure to system

109

parameters (muscle origin and insertion coordinates). Analytical sensitivity studies can
be constructed using the explicit form of the feedback gain solution provided in this
study.

Another area of future work is to simulate muscular neurological deficiencies and
their effect on the human motion. Muscle spasticity is considered to be one of the
common neurological deficiencies in which the muscle’s activation curves show
constant distribution (Bowsher et. al., 1992). There are two means to carry out such
research. The first approach involves the testing of subjects with pre-existing
neurological deficiencies. The results can then be used to establish general trends of
gain values associated with each deficiency. A second, more analytical approach,
involves the parameterization of the muscle’s activation curves in terms of a set of nodal
points. These nodal points can be altered to simulate different neurological deficiencies.
A numerical sensitivity study can then be developed to examine the effect of changing

nodal points on the estimated feedback gains.

Appendices

APPENDIX A

MUSCLE CONTRACTION DYNAMICS

The contraction dynamics of the muscle model pr0posed in Gordon et al. (1986) and
shown in Figure 2.2.2 is derived in this appendix. Since the volume of the muscle is
assumed to be constant, the muscle pennation angle at any muscle fiber length is related to

the pennation angle measured at muscle optimal length as follows
Lo- sinao = L,,,- sina = w = constatnt (A.1)
Differentiating Equation (A. l) with respect to time gives
a = (—L'm/L,,,) - tana (A2)

From geometry, the musculotendon length is related to the muscle fiber length at any

time as follows (see Figure 2.2.2.)

Lm, = L,+L,,,- cosa (A3)
where L, is the tendon length. After differentiating Equation (A3) and then combining it

with Equation (A.2), the musculotendon velocity V,,,, is related to the muscle fiber veloc-

ity as follow

V,,,, = V, + Vm/cosa (A.4)

where V, and V", are the tendon and muscle fiber velocities, respectively. Similarly
V", = V” + V“ (A.5)
where V“ and V“ are the series element and contractile element velocities, respectively.

By inspection of Figure 2.2.2, the force in the tendon p is related to the muscle fiber

110

111

force, pm as follows
p = pm - cos(0t) (A.6)

Note that p is the force passing through the muscle tendon. Since experimental data

of the tendon are given in terms of its stiffness K which is defined mathematically as

K = fif/V, (A.7)

the differentiation of Equation (A.6) with respect to time yields

dp _ dpm. - . - . -
E _ a; cos(0t) pm a. srn(0t) (A3)

Substitute from Equation (A2) for (x in the above equation yields

dP_"Pm. .5.
H7 — :1? cos(0t)+p,,, L

tan(0t) - sin(or) (A.9)

But from the equilibrium condition on the model, the following relation can be easily

verified

113... _ ff“ +515,»

dt _ dt dt (A.10)

Using the same stiffness to rate of force relation given in Equation (A.7), Equation

(A.10) can be expressed as

d
1’“ =K,,-V,,+Kp,-V

it (A.1 l)

pe
From the geometry of the muscle model shown in Figure (2.2.2), it is easy to note that

Lpe = Lm, hence the velocities are equal. Using the previous statement and Equation

(A.4), Equation (A.11) reduces to

 

112

{I}...

d, = (Ks, + K") - vm — K5,. V6, (A.12)

Substituting Equation (A.6) and Equation (A. 12) in Equation (A.9) yields

dp- . . - . . _P_. . 2
E — (K,,+K,,,) Vm cos(0t) K“ V“ cos(a)+(Lm) V”, (tan(0t)) (A.13)

By letting K ,,, = K S, + K p, and using Equations (AA and A.7), and with some alge-

braic manipulations, Equation (A.13) takes the form

dp_ K~Kma-cosa { K” V }
CC

 

:17 - K+K,,,O,- cosa. ”Fm. (A'14)
where
p 2
Km, = Km ~ cos(a) +i— - (tan((1)) (A.15)
m
In Equation (A.14), the musculotendon K n, K pa, cos(a) , L,,,, and V“ can be

evaluated using the constitutive equations given in Section 2.2.2. All terms are then

expressed in terms of Lm, , p , and a(t) , as in the case of V by using the geometric

CC ’
relations in this appendix. Thus the general form of Equation (A. 14), which represents

the muscle contraction dynamics, can be expressed as

f}? = m... V....p.a(t)) (A16)

APPENDIX B

BLOCK OBSREVER
STATE SPACE REALIZATION

Kailath (1980) gives a detailed account of various state space realization of multi-
input multi-output (MIMO) systems. In this appendix we shall introduce only one of the
MIMO state spaces, the Block Observer form of an Autoregressive (ARX) input/output
discreet model.

Consider the following general second order ARX input/output model. This equiva-

lent to n = 2 in Equations (2.1.5; 2.3.3),
(MEI ' 9" +712 - 9’2) - C(q) = (52 ' 44“? - 9" +52) ° W) 03-1)
where C( q) is the output vector, C(q) e R' , §(q) is the input vector, §(q) e Rm , and

Xi: A; 3,, and, B; are constant coefficient matrices. It can easily be shown that the

Block Observer state space representation of the system in Equation (B. 1) takes the form

(Strejc, 1981)

X(k+1) = A-X(k)+B-§(k) (B.2)
C(k) = C-X(k)+D'§(k) (3.3)
where
A - ’5 I
-A2 [0] 2rx2r

113

and

114

—*
Bl
—*
2 2er

APPENDIX C

EQUATIONS OF MOTION

The dynamic equations of motion of the 4 DOF skeletal model shown in Figure 3.1.2

are given in a vector matrix form in Equation (3.1.9) of Section 3.1

 

 

 

 

J(e)-é =B(e)-62+C(e)+D-M(e,).P (3.1.9)
Where
all —a12cos(01 - 02) a13cos(0, - 03) -a14cos(01- 04)
1(9) _ -a12cos(0l - 02) a22 -a23cos(93 — 62) a24cos(64 — 02)
a,3,cos(91 - 63) —a23cos(03 - 02) a33 —a34cos(03 — 04)
L—a14cos(01 — 64) aucos(04 - 02) ~034cos(03 - 04) a“ _
(C.l)
' o —alzsin(91—02) -a,3sin(61-03) -al4sin(61—04)
8(9) _ -alzsin(91- 92) 0 —a23 sin(93 — 02) ausin(02 - 04)
.0143in(01-64) —0248in(62-94) "d348in(e3 -64) 0 4
(C2)
p-clgcosfl;
C(6) = czgmsez ((2.3)
-c3gcose3
c4gc0804

 

 

115

 

116

F200-
D —2—20
022
_00—2

.1

(CA)

 

 

and

2 2

“12 = 62°11
“13 = ca'lr
“14 = 04'11

2 2
“23 =03'l2
024 = C4'12
2 2
a33=I3+M3'lc3 +m4‘l 3

“34 = 04°13

2
a44=I4+m4-lc4

c1: ml-lcl+(m2+m3+m4)~ll
02 = mz-lc2+(m3+m4)-lz
C3 = m3°lc3+m4’l3
C4 =m4’lc4

The parameters given in the above equations are defined as follows

 

117

I, = moment of inertia of segment i about its center of mass.

1, = length of segment i.

la = distance of body segment center of mass from distal end.
mi: mass of segment i.

9i = the angle that segment i makes with the horizontal.

6a, 01,, and 9h are the ankle, knee, and hip joint angles, respectively.
91T=[ 0,, 0k, 0h]T; joint angles vector.

6T: [ 01, 02, 03, 04]T; system output vector.

The linearization of Equation (3.1.9) is tacking place around the operating points
(nominal data) obtained and used in the identification of the closed loop system. The
method of linearization defined in Section 2.2.3 is dependent on the evaluation of the Jaco-
bian matrices in terms of the system states (positions, velocities) and system inputs. Let

us rewrite Equation (3.1.9) in the following form

9 = f(& 9.1”) (C5)

where
no, 9, P) = J"‘(e).(19((~))-e2 + C(e) + D . M(e,) - P) (C6)

The function f ( 6, O, P) e R4 , is a nonlinear function of system states. The gradient

of f with respect to 0,- e 9 is defined as the differentiation of the function with respect to

G.

 

118

-1
V9, f =% (e) - (3(9) . 92+C(9)+D-M(6,) . P)

_1 B - 2 3C 3M
+ J (9) - (SEQ) - 9 4.397(9) +D . 39—561). P) (C.7)

1

Note that if I -1 (6) exists then 3%: (6) can be computed from the following for-
i

mula

aJ"

_ -1 _BJ . —1
36‘, (e) _ J (9) RIG) 1 (9) (C3)

-1
Thus the computational difficulty of finding 3‘;— (O) is avoided by using Equation
1'

(CS). On the right hand side of Equation (C7), the moment arm matrix is not explicitly

defined in terms of segmental angles but is computed by using inter-segmental angles,

joint angles, as shown in Section 2.2.1. Hence, the term gig—{(61) is computed from
i

3&9» = 3’59.) --:—g‘-:+§%(e,) --:-69-’:+§—9‘4’;(9,) g3: (C9)
where
9,, = (62—6,)—1t/2 (C10)
9,, = 92—9, (C11)
9;. = 94-93 ((2.12)

as, 39,, as,

The terms 30,330? midi-9: 1n Equation (C.9) takes one of the values (1, 0, -1)

119

depending on which segmental angle is used in the differentiation. The Jacobian matrix

evaluated with respect to the vector 8 is then expressed as
Ve f = [V9, f V92 f V0, f V9, f] (C.13)

Similarly, the Jacobian matrix of f with respect to 0,- is evaluated from Equation (C.6)

as follows

. 2
V. f = J“(e) - 8(6) 3% ((2.14)

I

Therefore the whole J acobian matrix V9 f reduces to

V,3 f = 2 - J“(9) - 3(9) - diag(0,-) (C.15)

Finally, Vp f matrix is evaluated following the same procedure used to obtain Equa-

tion (C.15) and can be expressed as follows (C.16)
V,» f = J"(9) - D - M(9,) . 17 (C17)

and I7 is a (7x 7) identity matrix.

 

References

 

10

11

12

13

REFERENCES

Abdel-Aziz Y.I. and Karara H.M.: Direct Linear Transformation from Comparator
Coordinates into Object Space Coordinates in Close-Range Photogrammetry, Pro-
ceedings of the Symposium on Close-Range Photogrammetry, January 26-29, 1971.

Alexander R.McN., and Vernon A.: The Dimensions of Knee and Ankle Muscles and
the Force they Exert, J. Human Movement Studies, 1, pp 115, 1975.

Altringham J .D. and Johunson I.A.: The pCa-tension and velocity characteristics of

skinned fibers isolated from fish fast and slow muscle, J. Physiol.(Lond), 333, p 42,
1982.

Agarwal G.C., Berrnan B.M. and Stark L.: Studies in Postural Control Systems 1.
Torque disturbance input, IEEE Trans. Syst. Soci. Cybem. $506, 2, pp 116, 1970.

Anderson III B.A.: Storage and Utilization of Elastic Strain Energy During Human
Jumping: An Analysis based on Predictions of an Experimentally Verified Optimal
Control Model, Masters Thesis, University of Texas at Austin, 1992.

Bahler A.S.: Series Elastic Component of Mammalian Skeletal Muscle, Amer., J.
Physiol., 213. PP 1550, 1967.

Bahler A.S., Fales, LT. and Zieler K.L.: Dynamic Properties of Mammalian Skeletal
Muscle, Gen. Physiol, 51, pp369, 1968.

Barker R.A.: Neuroscience an illustrated guide, Ellis Horwood, N.Y., 1991.

Baumann W.T. and Rugh W.J.: Feedback Control of Nonlinear Systems by Extended
Linearization, IEEE Trans. on Automatic Control AC-3I, 1, pp 40, January 1986.

Becket R. and Chang K.: An evaluation of the Kinematics of Gait by Minimum
Energy, J. Biomechanics, 1, pp 147, 1968.

Belanger A.Y., McComas AJ. and Elder G.B.C: Physiological Properties of Two
Antagonistic Human Muscle Groups, Eur: TI Appl., Physiol. (51), p. 381, 1983.

Bemis T.A. : The parallel squat exercise: A Biomechanical comparison between nor-
mal and anterior cruciate ligament deficient knees, Ph.D. Thesis, Department of
Biomechanics, Michigan State University, 1992.

Bernstein N.A.: On The Motion Synthesis, Medgiz, Moscow, 1947.

120

 

14

15

l6

17

18

19

20

21

22

23

24

25

26

121

Bowsher K.A., Damiano D.L. and Vaughan C.L.: Joint Torques and Co-Contarctions
During Gait for Normal and Cereberal Palsy Children, Proceeding of NACOB II,
The Second North American Congress on Biomechanics, Chicago, August 24th-
28th, 1992.

Boyed I.A.: The isolated mammalian muscle spindle, Trends Nuerosci., 3, pp 258,
1980.

Boyed LA. and Roberts T.D.M.: Proprioceptive discharges from stretch receptors in
knee joint of the cat, J. Physil., 122, pp 38, 1953.

Brain K.: Evaluation of Generalized Model of Human Postural Dynamics and Con-
trol in Sagittal Plane, Biol. Cybern., 61, pp 37, 1989.

Brand R.A., Crowninshield R.D., Wittstock C.E., Pedersen D.R., Clark CR. and Van
Kricken F.M.: A Model of Lower Extremity Muscular Anatomy, J. Biomed. Eng,
104. PP 304, 1982.

Brooks V.B.: The Neural Basis of Motor Control, Oxford University Press, NY.
1986.

Bruker R.E., Levine D.N., Tsains P. and Zajac R.E.: Physiological types and his-
tochemical profiles in motor unit of cat gastrocnemius, J. Physiol., 234, p 273, 1973.

Camana P.C., Hemami, H. and Stockwell C.W.: Determination of Feedback for
human: An Optimal Control Model for Analyzing Human Postural Balance, J.
Cybem., 7, pp 199, 1977.

Carlson ED: Kinematic Studies on the Mechanical Properties of Muscle, In Tissue
Elasticity, Remington, J.W. Washington, Amer Physiol. Soc., 1957.

Chance M.A. and Bayazitolque V.O.: Development and Application of Generalized
d-Embark Force for Multi-Freedom Mechanical Systems, ASME Trans, J. Eng. Ind.,
93. PP 317, 1971.

Chao E. Y. and Rim K.: Application of Optimization Principle in Determining the
Applied Moments, J. Biomechanics, 6, pp 497, 1973.

Chen C.T.: Introduction to Linear Control Theory, Holt Rinehart and Winston, N.Y.,
1970.

Chow GK. and Jacobson D.H.: Studies of Human Locomotion Vra Optimal Pro-
grarning, Math. Biosci., 10, pp 239, 1971.

 

27

28

29

30

31

32

33

35

36

37

38

39

122

Clauser C.E., McConville J.T. and Young J.W.: Weight, Volume, and Center of mass
of segments of human body, Aerospace Medical Research Laboratory, AMRL-TR-
69-70, 1969.

Close R.I.: Dynamic Properties of Mammalian Skeletal Muscle, Physiol. Rev., 52,
pp129,1972.

Collins E.G., Phillips DJ. and Hyland D.C.: Robust Decentralized Control Laws for
ACES Structure, IEEE Control System Magazine, pp. 62, April 1991.

Crago P.E., Chizeck H.J., and Hambrecht F.T.: Sensors for use with functional neu-
romuscular stimulation, IEEE Trans. Biomed. Eng, Vol. 33, p. 256, 1986.

Crowe A.: Mechanical Model of Muscle and its Application to the intrafusal Fibers
of the mammalian Muscle Spindle, J. Biomechanics, 3, pp 538, 1970.

Crowninshield R.D.: Use of Optimization Techniques to Predict Muscle Forces,
ASME Trans. J. Biomed. Eng., 100, pp 88, 1978.

Crowninshield RD. and Brand R.A.: A Physiological Based Criterion of muscle
force Prediction in locomotion, J. Biomechanics, 14, 11, pp 793, 1981.

Dempster W.T.: Space Requirements of the Seated Operator, WADC Tech. Rept. p.
55, Wright. Patterson Air Force Base, Ohio, 1955.

Dhaher Y.Y.: Characterization of the position feedback synergy in uni- and bi- artic-
ular mucsle in the skeletomotor, First Annual Conference of The International F ES
Society (IFESS ' 96), Cleveland FES Center, Cleveland OH, May 14-16, 1996.

Dhaher Y.Y. and Soutas-Little R.W.: Postural Balance Modeled as a Double Inverted
Pendulum, First Annual North American Clinical Gait Laboratory Conf., Portland,
Oregon, April (7-9), 1994.

Diener HQ and Dichgans 1.: On the role of Vestibular, Visual and Somatosensory
information for Dynamic Postural Control, Progress in Brain Research, (Pompeiano
O. and Allum J.H.J., Ed), Elsevier, Amesterdam, 76, pp 253, 1988.

Ebashi S. and Endo M.: Calcium Ion and Muscle Contraction, Prog. Biophy, 18,
pp.l23,1968.

Ettrnema G.J.C. and Huijing P.A.: Architecture and elastic properties of the series
elastic element of muscle - tendon complex, In Multiple Muscle Systems: Biome-
chanics and Movement Organization, (Winters J.M. and Woo S.L-Y., Ed), Springer-
Verlage, New York, 1990.

4o

41

42

43

45

46

47

48

49

50

51

52

53

123

Fallside E: Control System Design by Pole-Zero Assignment, Academic Press, Sep-
tember 1974.

Finley ER. and Karpovich P.V.: Electrogoniometric analysis of normal and patho-
logical gaits, Res-Quart, 35, p 379, 1964.

Fisher, 0.: Theoretische Grundlayen Fureine Mechanik der Lebenden Korper, pp 52,
Teubner, Berlin, 1906.

Franklen H.M., Veltink P.H., Baardman G., Redmeijer R.A. and Boom A.B.K.:
Cycle-to-Cycle control of the swing phase of paraplegic gait induced by functional
electrical stimulation, Med. Biol. Eng. Comput. in press.

Franklen H.M., Veltink P.H., Trjsmans R., Nijmeijer H. and Boom A.B.K.: Identifi-
cation of the Quadriceps-Shank dynamics using Randomized Interpulse internal
Simulation, IEEE Trans. Rehab. Eng, 3, 2, June 1995.

Fuglevand A.J., Winter, D.A., Patla A.E., and Stashuk D.: Detection of Motor Unit
action potentials with Surface electrodes: influence of electrode size and spacing,
Biol. Cybem. (67). P. 143, 1992.

Fung Y.C.B.: Elasticity of Soft Tissues in Simple Elongation, Am. J. Physiol., 213,
pp 1532, 1967.

Gill RE. and Murray W.H.: Practical Optimization, Academic Press, London, 1981.

Gordon M.E.: An anlaysis of the bimechanics and muscular synergeies of human
standing, Ph.D. thesis, Department of Mechanical Engineering, Stanford University,
September 1990.

Goh CJ. and Teo K.L.: Control Parametrization: A unified approach to optimal con-
trol problems with general constraints, Automatica, 24, pp 3, 1988.

Goldstein H.: Classical Mechanics, Addison-Wesley, Mass., 1959.

Golliday G.L.Jr. and Hemami, H.: Postural Stability of the tow degree of freedom
Biped by general linear feedback, IEEE Trans. on Automatic Control, February, pp
74, 1976.

Gordon M.E., Hoy M.G., Zajac RE, and Maclean K.E.: The Musculoskeletal Model
of Human Lower Extremity, RESNA 9th Annual Conference, Minneapolis, Minne-
sota, pp 448, 1986.

Gordon A.M., Huxley A.F. and Julian F.I.: Variation in the Isometric Tension with
Sarcomere length in vertebrate muscle fibers, J. Physiol., 184, p 170, 1966.

 

54

55

56

57

58

59

61

62

63

65

67

68

124

Green D.G.: A Note on the Modeling Muscle Physiological Regulator, Med. Biol.
Eng., 7, pp 41, 1969.

Griffin J.E., Karselis TC. and Currier DP: Physical Agents for Physical Therapists,
Charles C Thomas Publisher, 1982.

Grigg P. and Greenspan B.J.: Response of Primate Joint Affemet Neurons to
Mechanical Stimulation of Knee Joint, J. Neurophysiol. 4, 1, 1977.

Hagedom P.: On The Stability of Steady Motions in Free and Restrained Dynamical
Systems, J. App. Mech., 46, pp 427, 1979.

Hanavan E.P.: A mathematical model of human body, Aerospace Medical Research
Laboratories, AMRL—TR-64-102, 1964.

Hasan Z.: A Model of Spindle Affemet Response to Muscle Stretch, J. Neurophysi-
ology, 49, 4, April 1983.

Hatze H.: Theory of Contraction and Mathematical Model of Striated Muscle, J.
Theory. Biol, 40, p 219, 1973.

Hatze H.: Biodynamic Models and Their Applications, J. Acoust. Soc. Am., 50, pp
1397, 1971.

Hatze H.: The Complete Optimization of Human Motion, Math. Biosci, 28, pp 99,
1976.

Hatze H.: A Control Model of Skeletal Muscle, Ph.D Thesis, University of South
Africa, January 1975.

Hatze H.: Zum Thema: Eine Mechanik der Leibesubungon, Leibesub. - Leibeserz.,
6, PP 6, 1965.

Hatze H.: Die Biophoronome Erfassung der Menschlichen Bewegung in Einem Sys-
tem Particular Differentialgleichunger, Leibesub. - Leibeserz., 7, pp 3, 1967.

Hatze H.: A myocybemetic control model of skeletal muscle, Biol. Cybern., 25, pp.
103, 1977.

Hatze H.: The Change - Transfer Model of Myofilaments Interaction: prediction of
Forces Enhancement and Related Myodynamic Phenomenon, In Multiple Muscle
Systems: Biomechanics and Movement Organization, (Mnters J .M. and Woo S.L-
Y., Ed), Springer-Verlage, New York, 1990.

He J ., Levine, and Loed G.E.: Feedback gains for correcting small perturbations to
statnding posture, IEEE Trans. on Automatic Control, Vol. 36, p. 322, 1991.

69

70

71

72

73

74

75

76

77

78

79

80

81

82

125

Helal J-N. and Bouissou P.: The Spacial Integration effect of Surface Electrode
Detecting Myoelectric Signal, IEEE Trans. on Biomed. Eng. (39), (11), p. 1161,
November 1992.

Hemami H.: Modeling, Control, and Stimulation of Human Movement, CRC Criti-
cal Reviews in Biomedical Engineering, 13, 1, 1985.

Hemami, H. and Camana P.C.: Nonlinear Feedback in Simple Locomotion Systems,
IEEE Trans. Automatic Control AC-ZI, 6, pp 855, 1976.

Hemami, H. and Cvetkovic V.S.: Postural Stability of Biped Models via Lyapunove’s
Second Method, IEEE Trans. Automatic Control A022, 1, pp 66, 1977.

Hemami, H. and Golliday G.L.Jr.: The Inverted Pendulum and Biped Stability, Math.
Biosc., 34, pp 95, 1977.

Herman R.N., Grillber S., Stein P.S.G. and Stuart D.S.: Neural Control of Locomo-
tion, Plenum Press, N.Y., 1976.

Hemami H. and Jaswa V.C.: On a Three - Link Model of the Dynamics of Standing
up and Setting down, IEEE, Trans. System, Man, Cybem, SMC-8, pp 115, 1978.

Hill A.V.: The Heat of Shortening and the Dynamic Constants of Muscle, Proc. Roy.
Soci. B. (Land), 126. pp 136, 1938.

Hill A.V.: The Series Elastic Element Component of Muscles, Proc. Roy. Soci. B.
(Loud),137, pp 273, 1950.

Hill A.V.: The Mechanics of Active Muscle, Proc. Roy. Soci. B. (Land), 141, pp 104,
1953.

Hill D.K.: Tension Due to Interaction Between the Sliding Filaments in Resting
Striad Muscle, J. Physiology, 199, pp 637, 1968.

Horak F.B., Nashner L.M. and Diener H.C.: Postural Strategies associated with the
Somatosensory and Vestibular Loss, Exp. Brain Res., 82, pp 167, 1990.

Houk J .C., Singer J .J . and Goldman M.R.: An evaluation of length and force feed-
back to soleus muscles of decerebrated cats, J. Neurophysiol, 33, pp 784, 1970.

Hout RC. and Little R.W.: A Constitutive Equation for Collagen Fibers, J. Biome-
chanics, 5, pp423, 1972.

83

84

85

86

87

88

89

91

92

93

94

95

96

126

Hoy G.M., Zajac RE. and Gordon ME: A Musculoskeletal Model of human Lower
Extremity: The effect of muscle, tendon and moment arm on the moment-angle rela-
tionship of the Musclotendon actuators at the hip, Knee, and Ankle, J. Biomechan-
ics, 23, 2, p 157, 1990.

Huston R.L. and Passarello, CE: The Mechanics of Human Motion, In Human Body

Dynamics Impact, Occupational and Athletic Aspects (Ghista D.N., Ed), Oxford
University Press, London, pp203, 1982.

Inman V.T. Ralston H.J. and Todd R: Human Walking (Lieberman J .C., Ed), Wil-
iams and Mlkins, Baltimore, MD, 1981.

Iqbal K. and Hemami, H.: Stability and Control of frontal Four Link Biped System,
IEEE Trans. on Biomed. Eng., 40, 1, October 1993.

Isermann R.: Digital Control Systems Volume II, Springer-Verlag, 1989.

Ishida A. and Miyazaki 8.: Maximum Likelihood Identification of Postural Control
System, IEEE Trans. on Biomedical Engineering BME-34, 1, pp 1, January 1987.

Jenkins G.M., and Watts D.G.: Spectral Analysis and its Applications, San Francisco,
Ca, Holden - Day, 1969.

Jhonson G.A., Choi, N.Y., Tramaglini D.M., McMahon PJ. and Woo S. L-Y.: Vis-
coelastic Properties of Human Patellar Tendon: Age Related Changes, Proceeding of
NACOBII, the Zed North American Congress on Biomechanics, Chicago, August
24th - 28th, 1992.

Johansson R. and Magnussson M.: Human Postural Dynamics, CRC Crit. Rev. in
Biomed., 18, 6, 1991.

Kailath T.: Linear Systems, Prentice-Hall, Englewood Cliffs, N.J., 1980.

Kane TR. and Scher M.R.: Human Self - Rotation by Means of Limb Movement, J.
Biomechanics, 3, pp 39, 1970.

Khang G. and Zajac F.E.: Paraplegic Standing Controlled by Functional Neuromus-
cular Stimulation Part I - Computer Model and Control-System Design, IEEE Trans.
on Biomed. Eng., 36, 9, pp 873, September 1989.

Khang G. and Zajac F.E.: Paraplegic Standing Controlled by Functional Neuromus-
cular Stimulation Part II - Computer Model and Control-System Design, IEEE
Trans. on Biomed. Eng., 36, 9, pp 885, September 1989.

Koivo A.J. and Guo T.H.: Adaptive Linear Controller of Robotic Manipulators, IEEE
Trans. Automatic Control AC-28, 1, p 162, 1983.

 

97

98

99

100

101

102

103

104

105

106

107

108

109

110

127

Komistek R.D., Soutas - Little R.W., Stich J .B., Dhaher Y.Y., and Paxson R.D.:
Dynamic Model of Human Lower Extremity Joint Forces, Zed World Congress of
Biomechanics VRIJE UNIVERSITEIT, Amesterdam, The Netherlands, pp 10, July
1994.

Kralj A., Bajd T., and Thrk R.: Electrical stimulation providing functional use of
paraplegic patient muscle, Med. Prol. Technol., Vol. 7, p. 3, 1980.

Krogman WM. and Johnston F.E.: Human Mechanics, Four monographs a bridged.
Technology, Doc. Repot. No. AMRL - TOR - 63 - 123, Wright - Patterson Air Force
Base Ohio.

Kuo A.: An Optimal Control Model for Analyzing Human Postural Balance, IEEE
Trans. Biomed. Eng., 42, 1, pp 87, 1995.

Lawson CL. and Hanson R.J.: Solving Least Squares Problems, Prentice-Hall Inc.,
Englewood cliffs, New Jersey, 1974.

Lee C.S.G. and Chung M.L.: An adaptive control strategy for mechanical Manipula-
tors, IEEE Trans. Automatic Control AC—29, 9, p 837, 1984.

Lee C.S.G. and Chung M.L.: Adaptive perturbation control with feed forward com-
pensation for robot manipulators, Simulation, 44, 3, p 127, 1985.

Lehman S.L. and Stark L.W.: Perturbation analysis applied to eye, head, and arm
movement models, IEEE Trans. Syst. Man Cybem. SMC-I3, 5, pp 972, 1983.

Leu M.C. and Hemati N.: Automatic Symbolic Derivation of Dynamic Equations of
Motion for Robotic Manipulators, ASME J. Dynamics, System, Meas., and Cont,
108, ppl72, Sept. 1986.

Lewis L.F.: Applied Optimal Control and Estimation: Digital design and estima-
tion, Prentice Hall, Englewood Cliffs, NJ, 1992.

Levine W.S. and Zajac F.E.: Unpublished analysis, 1984.

Lieh J ., Rodgers M.M., Tummarkota S. and Schrag D.R.: Multibody Dynamics for
Biomechanical Systems, ASME Advances in Bioengineering BED-22, pp 597, 1992.

Ljung L.: System Identification - Theory for Users, Egelewood Cliffs, NJ, Prentice-
Hall, 1986.

Magnus R.: Some Results of Studies in the Physiology of Posture, Cameron Prize
Lecturers, Lancinate, 211, 2, pp 531, 1926.

 

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

128

Magnussson M.: On The OptoKinetic Mechanism, Ph.D. Thesis, Department of
Oto-rhinolaryngology, University of Lund, 1986.

Markhede G. and Grimby 6.: Measurements of strength of the hip joint muscles,
Scand. J. Rehab. Med. (12). P. 169, 1980.

Maros D. and Orlanda N .: Contribution to the determination of the equations of
motion for multi-degree of freedom Systems, ASME Trans, J. Eng. Ind., 93, pp 191,
1971.

Nagurka ML. and Yen V.: Fourier-Based Optimal Control of Nonlinear Dynamic
Systems, ASME Trans. J. Dyn. Syst., Meas. and Cont, 112, pp 17, 1990.

Nashner L.M.: A Model describing Vestibular detection of body Sway motion, Octa
Otolaryngology, 72, pp 429, 1971.

Nashner L.M.: Vestibular Postural Control Model, Kybemetik, 10, pp 106, 1972.

Nashner L.M.: Vestibular and reflex Control of Normal Standing, Advances in
Behavioral Biology, 7, (Stein R.B., Pearson K.G., Smith RS. and Redford J .B., Ed),
Plenum Press, N.Y., 1973.

Needham D.M.: Macina Camits, Cambridge University Press, 1971.

Nubar V. and Contini R.: A minimal principle in biomechanics, Bull. Math. Biophys.
23. pp 337, 1961.
Ortega J .M.: Matrix Theory - A second course, Plenum Press, 1987.

Owens D.H.: Multivariable and Optimal Systems, Academic Press, London, 1981.

Pandy M.G., Anderson EG. and Hull D.G.: A Parameter Optimization Approach for
the Optimal Control of A Large-Scale Musculoskeletal System, ASME Trans. J.
Biom. Eng., 114, pp 450, November 1992.

Pandy M.G., Garner BA. and Anderson F.C.: Optimal Control of Non-ballistic Mus-
cular movements: A Constraint-Based Performance Criterion for raising from Chair,
ASME Trans. J. Biom. Eng., 117, pp 15, February 1995.

Pandy M.G., Zajac F.E., Sim E. and Levine W.S.: An Optimal Control Model for
Maximum-High Human Jumping, J. Biomechanics, 23, 12, pp 1185, 1990.

Pierrynowski MR. and Morrison J .B.: A Physiological Model of the evaluation of
Muscular Forces in Human Locomotion: Theoretical Aspect, Math. Biosci, 75, pp
69, 1985.

 

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

129

Pringle J .W.S.: Models of Muscle, Symposium Soci. Exper: Biol, 1, pp41, 1960.

Raibert M.H. and Graig J .J .: Hybrid Position/Force Control of Manipulators, ASME
Trans. J. Dyn. Syst. Meas. and Control, 103, pp 126, 1981.

Reichardt W.: Analogy between noiogram information and computation of relative
movement by visual system of the fly, Naturwissenschaften, 67, pp 8411, 1980.

Reichardt W. and Poggio T.: Figure ground discrimination by relative movement in
visual system of the fly 1. Experimental Result, Biol. Cybem. 35, pp 81, 1979.

Rugh W.J.: Design of nonlinear compensator for nonlinear systems by extended lin-
earization technique, in Proc. IEEE Conf Decision Control, Las Vegas, pp 69, NV,
1984.

Sale D., Quinlan J ., Marsh E., McComas A.J. and Belanger A.Y.: Influence of Joint
position on ankle plantar flexion in humans, J. Appl. Physiol. Respirat. Environ.
Exer: Physiol. (52), p. 1636, 1982.

Seirge A. and Arvikar R.L.: A mathematical Model for the Evaluation of Forces in
Lower Extremities of the Musculoskeletal System, J. Biomechanics, 6, pp 313, 1973.

Seirge A. and Arvikar R.L.: Biomechanical Analysis of Musculoskeletal Structure
for Medicine and Sports, Hemisphere Publishing Corporation, New York, 1989.

Scudder G.N.: Torque curves produced at the knee during isometric and isokinetic
exercises, Arch. Physiol. Med. Rehab. (61), p. 68, 1980.

Shiavi R.: Electromyographic patterns in adult locomotion: A Comprehensive
Review, J. Rehabilitation Research and Development, 22, p 573, 1985.

Shiping Ma. and Zahalak 6.1.: Activation dynamics for a distribution-moment
model, Math. Comput. Modelling, 11, pp 778, 1988.

Sirisena HR. and Tan K.S.: Computation of Constrained Optimal Controls Using
Parametric Techniques, IEEE Trans. J. Automatic Control AC-19, pp 431, 1974.

Soderstrom T.: Lecture Notes in Identification, Control System Analysis Group,
Uppsala University, Sweden, 1984.

Soderstrom T. and Stoica P.: System identification, Prentice-Hall International
Series in Systems and Control Engineering, 1989.

Sonnenblick E.H.: Series Elastic and Contractile elements in Heart Muscle Changes
in Length, Annual Meeting of the American Physiol. Socie., 14th April 1964.

 

141

142

143

144

145

146

147

148

149

150

151

152

153

130

Soong T.T. and Huang W.N.: A Stochastic Model for Biological Tissue Elasticity in
Simple elongation, J. Biomechanics, 6, pp451, 1973.

Stanic U. and kaoczy A.: Closed loop positioning of herniplegic patient’s joint by
means of functional electrical stimulation, IEEE Trans. Biomed. Eng. Vol. BME-ZI,
p.365,1974.

Strejc V.: State Space Theory of Discrete Linear Control, Prague, Czechoslovak,
John Wiley & Sons, 1981.

Sugi H.: Tension Changes During and After Stretch in frog Muscle Fiber, J. Physiol.,
134, Pp 237, 1972.

Takeyaki M. and Arimoto S.A.: A new feedback method for dynamic control of
manipulators, ASME Trans. J. Dyn. Syst. Meas. and Control, 102, pp 119, 1981.

Tashman S., Zajac FE, and Perkash 1.: Modeling and Stimulation of Paraplegic
Ambulation in Reciprocating Gait Orthosis, ASME Trans. J. Biomed. Eng., 117, pp
300, August 1995. '

 

Tewell BR. and Wilkie D.R.: An Analysis of the Mechanical Components of frog’s
Striad Muscle, J. Physiol., 143, pp 515, 1958.

Vane J .M. and Sitchin A.: Derivation of the first -order difference equations for
dynamical systems by direct application of hamiltonian principle, ASME Trans.,
J.Appl. Mech., 37, pp 276, 1970.

Von Gierke, H.E.: Biodynamic Models and Their Applications, J. Acoust. Soc. Am.,
50, pp 1397, 1971.

Vukobratrovic M., Borovac M., Surla D. and Stokic D.: Scientific Fundamentals of
Robotics 7: Biped Locomotion Dynamic, Stability, Control and Application,
Springer-Verlag, Berlin Heidelberg, 1990.

White S.C., Yack H.J., and Winter D.A.: A Three Dimensional Musculoskeletal
Model for Gait Analysis Anatomical Variability Estimate, J. Biomechanics, 22, pp
885, 1989.

Walters R.L., Frazier J., Garland D.E., Jordan C., Perry 1.: Electromyographic Gait
Analysis befor and after Operative Treatment for Henriplegic Equinus and Equino-
varus Deformity, J. Bone and Joint Surgery, 64-A, 2, pp 284, 1982.

Whitney D.R.: Force Feedback Control of manipulators fine movement, ASME
Trans. J. Dyn. Syst. Meas. and Control, 99, pp 91, 1977.

154

155

156

157

158

159

160

161

162

163

164

165

131

Wilkie D.R.: Measurement of the Series Elastic Component at Various times during
a Single Muscle Twitch, J. Physiol. (Lond), 134, pp 527, 1956.

Williams W.J.: A System-Oriented evaluation of the role of Joint receptors and other
afferents in position and motion sensory, CRC Crit. Rev. in Biomed. Eng., 7, 3,
1991.

VVrnter D.: Biomechanics and Motor Control of Human Movements, John Wiley
&Sons,1nc., 1990.

“Enter D., Rau G., Kadefors R., Broman H., and Deluca C.: Units, Terms, and Stan-
dards in the reporting EMG research, Intern. Soci. of Electrophysiology and Kinesi-
ology, August, 1980.

“Enters J .M. and Stark L.: Muscle Models: What is Gained and What is Lost by
Varying Model Complexity, Biol. Cybern., 55, pp 403, 1987.

Winters J .M., Fuglevand A.J., and Archer S.: Crosstalk in Surface electromyography
in agonist Muscles, 8th annual East Coast Clinical Gait Laboratories Conference,
May 5-8, Mayo Clinic, Rochester, MN 1993.

 

Wittenburg J .: Dynamics of Systems of Rigid Bodies, B.G. Teubner, Stuttgart, 1977.

Yamaguchi G.T., Sawa A.G.W., Moran D.W., Fessler M.J. and Winters J .M.: A sur-
vey of human Musclotendon actuator parameters, In Multiple Muscle Systems: Bio-
mechanics and Movement Organization, (Winters J .M. and Woo S.L-Y., Ed),
Springer-Verlage, New York, 1990.

Young L.R.: A Control Model of Vestibular System, International Federation of
Automatic Control (IFAC), Technical and Biological Problems in Control, A Cyber-
netic View, (Iberall S. and Reswick J .B., Ed), pp 543, Instrument Society of America,
Pittsburgh, 1970.

Zajac F.E.: Muscle and Tendon: Properties, Models, Scaling, and Application to Bio-
mechanics and Motor Control, CRC Critical Reviews in Biomedical Engineering,
17, 4. PP 349, 1989.

Zajac F.E., Topp EL. and Stevenson P.J.: Musclotendon Actuator Models for Use in
Compute Studies and Design of Neuromuscular Stimulation Systems, RESNA 9th
Annual Conference, Minneapolis, Minnesota, 1986.

Zhang Li-Q., Shiavi R. and \Vrlkes M.: Modeling and identification of human mus-
culoskeletal walking system, in Proceeding of the 22ed Southeastern Symposium on
System Theory, Cookeville, TN, pp 150, March 11 - 13, 1990.

 

132

166 Zhang Y-F. and Hemami, H.: Impact effects of Biped Contact with environment,
IEEE Trans. Syst. Man Cybem. SMC-l4, 3, pp 437, 1984.

167 Ziegler J .M., Anderson F.C., Pandy MG. and Whalen R.I.: Musculoskeletal Models
Computational Algorithms Supercomputers, ASME Advances in Biomedical Engi-
neering BED-22, pp 601, 1992.

168 Zipp B: Effect of electrode parameters on the band width of the surface e.m. g. Power
Spectrum, Med. Biol. Eng. Comput., 16, p. 537, 1978.

 

 

..L'l

"I71111111111111111111155