I'
1
A METHOD OF SYSTEM IDENTIFICATION APPLIED TO A
I
SIMULATION OF HUMAN EXERCISE
;i
by
Hugh Austin King
j B.S., Iowa State University, 1968
 M.S., New York University, 1970
[ M.D., University of Pennsylvania, 1977
A thesis submitted to the
f Faculty of the Graduate School of the
I
University of Colorado at Denver
1 in partial fulfillment
il
I of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
i 1992
I A* 3
I
This thesis for the Doctor of Philosophy
degree by
Hugh Austin King
has been approved for the
Department of
Mathematics
by
David Fox
Date
James Koehler
King, Hugh Austin (Ph.D., Applied Mathematics)
A Method of System Identification Applied to a Simulation of Human Exercise
Thesis directed by Associate Professor William Briggs
A simulation model of the inputoutput system consisting of a human riding
an exercise Hi cycle is developed. The input consists of work load and the out
put is oxygen consumption. The model is based on physiologic principles and
data collected on human subjects. Properties of the model are discussed and its
behavior is compared to that of human subjects under a variety of test inputs.
A method of, system identification is then developed. This method uses least
squares and its generalizations to estimate parameters in N ARM AX and Volterra
series models. A useful set of input test signals is determined and a method for
matching the properties of these signals to the system under test is discussed.
The identification method is first applied to a Wiener system consisting of a first
order linear subsystem cascaded with a clipper nonlinearity and the resulting
N ARM AX and Volterra series models are discussed. The identification method
is then applied to the simulation model and N ARM AX and Volterra models are
constructed. These models are compared and their agreement with the original
simulation is evaluated.
The form and content of this abstract are approved. I recommend its publication.
Signed
illiam Briggs
CONTENTS
i
J
INTRODUCTION
1 THE MODEL
The System, to Be Modeled
I
The Computer............
The Bicycle,............
The Measurement Devices .
,i
The Subject.............
The Muscle Subsystem
I
The Venous Subsystem
The Lung Subsystem .
The Cardiac Subsystem
The Visceral Subsystem
i
The Arterial Tree . .
The Overall System . .
An Example .............
1
6
6
6
8
13
16
18
34
38
39
41
41
42
43
Conclusion
50
53
53
53
59
62
73
73
80
81
81
83
93
94
98
100
103
130
136
139
142
142
ANALYSIS OF THE MODEL
{
Introduction ...........'..............................
Ramp Inputs............................................
Step Inputs............................................
Sine Wave Input........................................
Constant Switching Pace Inputs ........................
The Model with Noise...................................
Summary ...............................................
DEVELOPMENT OF THE IDENTIFICATION METHOD
Introduction ..........................................
The Data Record........................................
Candidate Models.......................................
The Volterra Series ...............................
The NARMAX Model...................................
The Assessment Rule ...................................
An Example ............................................
NARMAX Applied to the Example......................
Volterra Series Applied to the Example.............
Conclusion.............................................
IDENTIFICATION OF THE MODEL
Introduction
The NARMAX Model............................................146
i
The Volterra Model .................................172
Building Models in the Presence of Noise ...................180
5 SUMMARY AND FUTURE WORK 191
Summary of the Thesis.............................191
Future Work.......................................193
l
APPENDIX A 196
Wiener vs Hammerstein Models.....................196
BIBLIOGRAPHY
222
vi
ACKNOWLEDGEMENTS
The author would like to thank Carol Welch of the Denver Veterans Administra
tion Hospital Pulmonary Function Laboratory and Charles Irvin of the National
'I
Jewish Hospital Pulmonary Function Laboratory who allowed us to use exercise
l
testing equipment. Connie Borst and Paul Hogan, who work in the Veterans Ad
ministration Hospital Lab, and Steve Senn, Martin Carlos and Dave Gurka, who
i
work in the National Jewish Hospital Lab have been very helpful. The author
would also like to thank George Swanson, formerly of the University of Colorado
Health Sciences Center who also allowed us use of exercise testing equipment.
The author thanks his thesis committee and acknowledges the University of
Colorado at Denver College of Liberal Arts and Sciences, which supported some
of this work through a research assistantship.
INTRODUCTION
The determination of V02maxthe maximum rate at which a subject can
consume oxygenby graded exercise protocols has been in use for subjects rang
ing from elite athletes to victims of severe cardiopulmonary disability for a long
time [24]. In these graded protocols, the subject is presented with a gradually
increasing work load and the rate of oxygen consumption is monitored until ei
ther fatigue sets in, the subject cant maintain the load or the rate of oxygen
consumption plateaus. This maximum rate of oxygen consumption (V02max)
is then used both as a measure of conditioning state and to prescribe exercise
regimens.
/
In recent years two possibilities have led to further interest in this system.
One is the possibility of determining VOimax by submaximal tests. This is par
ticularly important in the post myocardial infarction patient or the patient with
severe cardiopulmonary limitations where one would prefer to avoid maximal
I
exercise. It is also important in limited, unusual environmentssuch as a space
I
craft before reentrywhere assessment of conditioning is important but, exercise
to maximal levels is inconvenient. The second possibility is that further insight
into both the conditioning state and the underlying physiologic system may be
obtained by studying the dynamic's of the system under inputs other than gradual
i
increase of the work load.
!
In this thesis we will consider a broader problem which includes all of the
abovethat off system identification for this input(work load)output(oxygen
consumption) system. The specific system that we will consider is a human riding
;!
an exercise bicycle. A similar analysis would apply to a treadmill system or an
il
arm cycle or any number of braked ergometer systems that exercise large muscle
:l
groups. We will use the bicycle system because there is a large installed base
I
of test equipment available in hospitals and our analysis will be applicable not
only to esoteric, expensive equipment in physiology labs, but also to equipment
in hospital cardiopulmonary evaluation units. (There is also a large number
of treadmills available, but their dynamic response is generally slower than the
bicycles and, thus, they are not as useful.)
f
The specific problem we consider is that of determining, from a given input
time series of work load verses time and the corresponding output of VO2 verses
time, how one caln best construct a model that represents the system. We take the
view that encoded in this mapping from the space of possible inputs to the space
i
of possible outputs is all the information about the system that we can know.
.j
Once this system identification has been done, determination of VO^max and the
study of system [dynamics is a relatively easy subproblem. More specifically we
I
will limit ourselve to building models where the output of the system is expressed
2
as a polynomial whose terms are products of previous system inputs and outputs;
:l
socalled lagged polynomials. As an example let the input be the time series x(t)
i
and the corresponding output be y(t). A polynomial model of the type we are
l
considering would be
ii
I,
y(i) = ax(t 1) + by(t 2) + cy(t 1 )x(t 2).
Such models have a long history in system identification [22]. If we restrict the
polynomials to have factors involving the input only then the method becomes the
impulse response method [25] in the linear system case and the Volterra method
[39] in the nonlinear system case. Similarly, if we allow the lagged polynomials to
include both inputs and previous outputs the method is the ARMAX model for
AutoRegressive Moving Average with eXogenous inputs [11] in the linear system
case and the NA'RMAX model for Nonlinear Autoregressive Moving Average with
eXogenous inputs [27] in the nonlinear case. We will develop methods for deter
mining the lagged polynomial model from inputoutput data. This will involve
determining appropriate inputs with which to test the system and determining
i
algorithms to construct the model from the resulting inputoutput pair.
The approach we will use is a combination of analysis and numerical ex
)
periment applied not to the original system, but to a simulation model of the
system that we will construct working in MAT LAB.1 By using this technique
i
it might appear that the result will be only as good as the model. We offer three
'MATLAB is a trademark of The Math Works, Inc.
3
responses to thi's: The first is that the simulation model matches the response of
I
the original system quite well, as' we will show. The second is that the experi
mental design questions we will addresssampling time, input shape, etc.are
i!
the same whether using a model or the original system and so the techniques
developed will be applicable even if the parameters or structures in the original
system are not i exactly those of the simulation. The third is that even if the
simulation had no relation to any physiologic system it is still a complex, well
J
defined system and its identification is a fascinating problem in its own right that
provides insightj into both system behavior and identification methods.
i
On the other hand, there are reasons for using a simulation instead of the
i!
original system." By using the simulation model we have been able to perform a
large quantity of experimentation that would have been impossible with human
:i
subjects. We have also been able to conduct experiments such as impulse re
sponse and sustained high work loads that could not have been done with human
i;
subjects. By using the simulation model we have access to observablessuch
I'
as true transfer functions, venous pool sizes and instantaneous oxygen con
centrations, and othersthat would not be available in a human experiment.
I
Finally, the experiments are repeatable. While human performance varies from
j
experiment to experiment and day to day, the simulation produces consistent
results.
The thesis is;! divided into five chapters. In the first we will develop the sim
ulation model. In the second we will evaluate the simulation model and present
qualitative and quantitative observations about the simulation and relate them
to the original human system. In the third chapter we will develop methods of
i
system identification and apply them to a simple system. In Chapter Four the
methods of Chapter Three are used to identify the simulated system. Finally, in
Chapter Five we will sumarize the results and present areas for future research.
I I
I
5
CHAPTER 1
THE MODEL
!!
' The System to Be Modeled
In this chapter we will build the simulation model. The physiological system
i
that we are modeling is a human riding an exercise bicycle. The input to this
physiological system is a work load determined by the braking force at the pedals
V
of the exercise bike and the output is the oxygen consumption of the subject
determined by measurement at the mouth. The physiological system, however,
is only a part of the entire system we must model. The major system blocks
are the bicycle,, the subject and the measuring devices as shown in Figure 1.1.
The experimental protocol is entered into a computer which then executes the
experiment and,'controls data collection. We will include the computer in our
i!
analysis and so the full system is as shown in Figure 1.2.
ii
I
The Computer
i
The data flow begins when the investigator either enters a protocol into the com
puter or chooses an established protocol already stored. This protocol consists
of levels of work load and times to switch from one level to the next. In some
i
available systems the number of switching times and levels is limited only by
6
Figure 1.1: Basic block diagram of experimental setup.
Figure 1.2: Block diagram of full experimental setup.
7
I
computer memory and the time between switches is limited only by the com
\
puter cycle time. However, in some commercial units these choices are limited by
the software sold with the package. We will explore the effect of the number of
i
I'
switches and the time between switches on the quality of identification in Chap
I
ter Three. The computer also administers the experiment in the sense that when
'i
I
the operator presses the start button and the experiment begins, the computer
'j
clocks the digitally stored protocol to a D/A (digital to analog) converter which
sends an analog braking signal to the bicycle. It also begins to sample the ports
of an A/D (analog to digital) converter that interfaces to the output data lines
for gas flow, gas concentration, heart rate and other variables of interest. The
computer does the necessary numerical processing such as data smoothing and
the calculation of oxygen volume from the measured values of total gas flow and
percentage of oxygen in the total gas. It also provides information so the operator
can follow the progress of the experiment. And finally, the computer stores the
data values for further analysis and provides graphical output, either online or
i
in later offline analyses.
The Bicycle
The bicycle receives the analog braking signal from the D/A converter. It has
its own mechanical and electrical properties (for example, inertia) and acts as a
subsystem with a known transfer function. The output of the bike subsystem is
presented to the] subject subsystem, as shown in Figure 1.2, as a force against
i 8
which the subject must work. The subject maintains a near constant pedaling
.1
rate and so this force translates into a work load with units of watts. (We will
use the MKS system of units.)
I
i'
We will model the bicycle as a linear system with known time constant rb
I'
(the time required for the system to reach 1 1/e of its final value during a step
;l
I
input.) This must be determined by measurement and varies widely from bike
i
to bike. In a high quality lab bike, 77, may be as low as 1 second while in some
commercial systems it is as high as 10 seconds. The equation governing the bike
is assumed to be
ny + y = i
(1.1)
where y(t) is the load presented to the subject and i(t) is the voltage presented
j:
to the bike to determine braking force. Integrating (1.1) gives:
il
y(t) = 7/(t0)exp ((t t0)fn) + (1 In) [ i(*) exp ((t s)/Tb)ds. (1.2)
In the simulation model there is an internal update time (or time step) that
is chosen to be small relative to all time constants of the system being modeled.
The switching times in the protocol are required to be multiples of this internal
time. Thus, during any of these time steps the input to the bike is constantij
:j
during the jth time step. The bike output for t in the jth time step then becomes
Vi*) = y{tj) exP ((* ~ tj)/n) + ij{ 1 exp (i tj)ln) (13)
where tj is the smarting time of the jth time step.
9
Figure 1.3: Input (solid line) and output (dashed line) of the bicycle.
Figure 1.3 shows the output of the bicycle (dashed line) corresponding to an
input (solid line) switching between 25 watts and 75 watts in a manner called
i i
pseudorandom,., a commonly used test input for identification of linear systems.
The internal time step for this run is 0.1 seconds and the time constant of the
bike is 5 seconds. The dashed line in the figure represents the bike output and
il
!
follows the rapidthenslow change expected from a first order linear system.
: 1
There is another system component that should be discussed with the
bicyclethe antialiasing filter. Since the data processing methods developed
,i
i'
later will be digital, we are dealing with a sampled data system and thus the
data must be filtered to avoid aliasing. (See [21] for a discussion of this topic.)
One way to, see the importance of this antialiasing filter in the setting of system
10
identification is as follows: Suppose we input a sine wave whose frequency is below
the Nyquist frequency1 into an arbitrary system and observe its effect. We have,
j
in a sense, identified the system at that frequency. Now add the first harmonic of
the sine wave to the input and suppose that it lies above the Nyquist frequency.
Without filtering the system will produce its response to the combined signal.
I
However, when ;the data are sampled, the energy at the higher frequency will
be folded onto a; lower frequency and the identification algorithm will mistakenly
il
interpret the resulting output as having come from the single sine wave below the
Nyquist frequency. Clearly this is an unsatisfactory result and so we must filter

the input as a function of the sampling frequency that will be used in the final
i
identification algorithm. On the other hand, it is known [29] that an input signal
that is used in an identification experiment must have enough frequency com
ponents and so we must be careful that we do not remove too many frequencies
with our filter. ;! We implement this filter in the simulation as a fourth order
Butterworth filter [42] placed early in the data path. Its effect is illustrated in
i
Figure 1.4 which ishows the effect of a varying sampling interval on an input that
switches every ten seconds between 25 watts and 75 watts. As one goes down
the figure the sampling time goes from 0.5 second to 1 second to 2 seconds and
finally to 5 seconds. One can literally see the high frequencies disappear from the
signal. Although the effect on the output of the system depends on the frequency
1The Nyquist frequency is the lowest frequency such that there are two samples per cycle
and, thus, its value in Hertz is one half the sampling rate in samples per second.
11
100
oli1''1
0 50 100 150 200 250
1 TIME(SECS)
Figure 1.4: Effect of sampling and filtering on an input. The input is a step
shown in the top plate. In the successive plates the sampling time increases and
the corresponding antialiasing filter becomes more narrow.
12
response of the components, it is clear, as we will discuss later, that this effect
is important. The actual realization of the filter in the true system can be in
i'
hardware between the computer and the bicycle or it can be realized digitally in
the computer. Causality is not required in these experiments because inputs are
i
!'
determined entirely before the experiment begins. The way we will implement
this filter in the simulation is with the wellknown signal processing method of
I
passing the input through the Butterworth filter and then reversing the output
I
and passing it again through the filter. (This is done for ease of implementation.
If we have a hardware filter or any filter with nonzero delay we could easily
correct for it.)
,1
i
The Measurement Devices
Returning to the block diagram of Figure 1.2 we will discuss the Measurement
Devices block. ln a typical experiment the raw variables that are measured
are instantaneous heart rate, respiratory rate, gas flow at the mouth, oxygen
consumption and carbon dioxide concentration at the mouth. Many systems
provide for simultaneous measurement of other physiologic observables such as
I
hemoglobin saturation by ear or finger oximetry.
The heart rate can be measured by counting the time between voltage peaks
of the QRS complex of the EKG and inverting this to give beat to beat heart
rate. In the systems we have used this value is available as a voltage output from
the EKG monitor and is input to the computer on one of the channels of the
13
A/D converter. .,The computer samples the channel and may or may not smooth
I
l(
it before storing' it depending on the experimental protocol.
ii
Gas flow is typically measured by a fan of low moment of inertia placed in
the gas flow path just distal to the mouth. The fan rotations are counted by a
photocell sensor that responds to reflected light pulses. The output current of the
'I
photocell is used to charge a capacitor that acts as an integrator. The capacitor
is discharged with each reversal of the fan and so the instantaneous voltage across
it corresponds to the volume of gas passing the fan so far in each breath. This
voltage across tile capacitor is also input to a channel on the A/D converter and
the computer samples it according to the experimental protocol. By measuring
il
the time between reversals, the computer can determine instantaneous respiratory
rate and by differencing this total volume signal the computer can determine the
volume of gas flow per sampling interval.
I
Some lab systems use a mass spectrometer to determine the composition of
the expired gas. These systems output analog signals corresponding to percentage
oxygen, percentage carbon dioxide and percentage nitrogen to channels on the
A/D converter. The percentage oxygen signal is sampled and multiplied by the
corresponding gals flow signal in the computer to yield oxygen flow. This signal
is integrated over an exhalation and compared to the value for the corresponding
inhalation to get the breathtobreath oxygen consumption. A similar calculation
l
may be done with the carbon dioxide data to get the breathtobreath carbon
dioxide production.
14
,1
Other systems perform carbon dioxide analysis with an infrared absorption
analyzer and oxygen analysis with a fuel cell. It is important to calibrate these
!:
measurement systems and to determine their inherent time delays. However, once
a
this is done we believe that they accurately reflect the physiologic variables and
we will model them in the simulation as a straight passthrough from the state
variables output, by the Subject Box in Figure 1.2.
i
We must mention, however, that the output from the system is selfsampled
on a breathtobreath basis. Although we have access to flow data, gas concentra
tion and other variables at arbitrarily high sampling rates (some systems sample
these in the kiloliertz range), these data are used only to calculate how much oxy
gen is inhaled arid how much is exhaled as discussed above. These are subtracted
to get one number per breath that represents the oxygen consumption during
that breath. We,have found that a subject can pace respirations to a nearly fixed
rate by varying only the depth of respiration within limits. There has also been
1
work done regarding best breathtobreath interpolation techniques. There are
many implementations of these methods and a real bike experimental system can
provide output time series sampled on a breath to breath basis or, by interpo
lation, at any sampling rate. Our simulation model will produce observations of
these variables as frequently as one observation per time step. I
I
I
15
The Subject
Now consider the Subject Box in Figure 1.2. The input to this box is the
work load presented by the braking force applied to the pedals and the output
is the vector of sampled physiologic observables as discussed in the previous sec
tion. Figure 1.5 shows the major anatomical components of the cardiorespiratory
system: L represents the lung and upper airways that are the source of external
oxygen for the system and the sink for carbon dioxide; C is the heart that pumps
1
the blood that carries the oxygen to the tissues; M represents the muscle com
partment; 0 represents the other consumers of oxygen (brain, viscera, etc.); V
represents the venous system that returns the oxygen depleted blood from the
muscles and viscera to the heart and lungs; and A is the arterial system that
carries the oxygenated blood from the heart and lungs to the peripheral com
partments. In this diagram one can trace the oxygen path from external air to
mitochondria that will eventually determine VO2 (the rate of oxygen consump
tion). U(t) is the work load input which will act on the muscles and the cardiac
subsystems as represented by the arrows in the figure. Y(t) is the observed out
put vector that acts as input to the Measuring Devices Box. We will now discuss
and model each subsystem in turn.
16
Y(t)
U(t)
Figure 1.5: Major anatomical components of the cardiorespiratory system.
17
Figure 1.6: Overall model of the muscle compartment.
ij
The Muscle Subsystem
For our purposes U(t) (exercise load) will affect the muscles in two ways. It
will cause oxygen to be locally consumed as glycolysis occurs to fuel the muscle
contractions. It will also regulate local blood flow to carry oxygen and other nu
j
trients into the muscle and to carry carbon dioxide, lactate and other metabolites
away. We will model these as two independent processesone with work load as
input and muscle blood flow as output and a second with work load as input
and muscle oxyjgen consumption as output as in Figure 1.6.
This type of model has been used before. ([5], [14].) A fundamental prop
j
erty of linear systems is superposition which implies that doubling input doubles
output. It is clear that the two processes discussed above are nonlinear for large
input variations, since if we repeatedly double the work load there will be a load
18
i
i
beyond which the flow no longer doubles. While the image is not so dramatic,
the local consumption of oxygen will also stop doubling as the load doubles be
yond some limit; It is also clear that these subsystems have dynamics in that the
response to an input is not instantaneousone does not note shortness of breath
i!
or increased heart rate immediately after increasing work load. Much research
has been directed toward determining the dynamics of these subsystems as dis
cussed in [5]. We will assume that they are first order. The simulation model is
modular and the implementation of higher order kinetics would be straight for
ward. One methpd for modeling nonlinear systems with dynamics is to isolate the
nonlinearity in one subsystem with no dynamics(ie. instantaneous response and
no memory) either preceded or followed by a second subsystem that models the
!i
dynamics and is(linear [22]. If the linear subsystem precedes the nonlinear static
subsystem the model is called a Wiener model and if the nonlinear subsystem
ij
precedes the linear subsystem the model is called a Hammerstein model as
shown in FigureJl.7.
We will represent both of the muscle level subsystems (work load to blood
flow and work load to oxygen consumption) in one or the other of these ways. In
the case of the blood flow subsystem, for the Wiener model the work load will be
input to a first order linear system. Since the simulation model is discrete in time
i
we will conceptually sample the input work load at each time step and the model
i
. l
will evolve between consecutive internal steps governed by the first order linear
model. The output of this subsystem will then become input for the nonlinearity
WIENER MODEL
HAMMERSTEIN MODEL
Figure 1.7: Wiener and Hammerstein models.
to yield the overall output. On a physiologic level, this means that the input work
load is delayed and its power density spectrum is modified by some underlying
mechanism. However, this is quite reasonable since we are seeing the modulation
of local vascular resistance presumably mediated by local messengers either
*!
chemical or neurologicalthat do not respond immediately or uniformly to work
f
load. This apparent force that is trying to alter flow is then modulated by the
reality that the arteriolar/capillary system has finite capacityas represented by
the horizontal asymptote on the function representing the nonlinear portion of
i
the model. This is seen in Figure 1.8 which shows such a Wiener system. In the
;l
Hammerstein model the cascade is reversed and the argument is that the input
work load is modulated by the nonlinearity first to yield an apparent force
ji
which becomes the input to the linear dynamic subsystem.
Similar heuristics apply to the oxygen consumption subsystem. In this case
the forces are not the opening of blood vessels that lead to changes in flow, but
20
Figure 1.8: Wiener model with limiting nonlinearity. LS is the linear system and
1
the saturating curve represents the nonlinear system.
are derived from chemical potentials that drive reactions down redox chains and
lead to changes jin oxygen consumption.
We must now decide which model to useWiener or Hammerstein. Figure 1.9
r
shows the result of a sine wave input with high frequency relative to the time
constants of the system to the Wiener model and Figure 1.10 shows the result
when the same sinusoid is input into the Hammerstein model. The nonlinearity
li
in each is the saturating type shown in Figure 1.8. The output in each case is
I,
1
nearly a sine wave with the same frequency as the input. However, the mean of the
output for the Hammerstein model is less than for the Wiener model. Figure 1.11
i
shows a similar plot for a sine wave input with low frequency relative to the time
constant for thel Wiener system and Figure 1.12 is the same for the Hammerstein
system. Note that here the outputs from both models are periodic and of similar
appearance, but' nonsinusoidal. Again the mean for the Hammerstein case is less
than for the Wiener case. Thus, although the models produce outputs that are
21
OUTPUT" UNITS
Figure 1.9: Output for high frequency sinusoidWiener model.
i
22
OUTPUT UNITS
I
0 50 100 150 200 250 300
TIME(SECS)
Figure 1.10: Output for high frequency sinusoidHammerstein model.
23
OUTPUT UNITS
46
0 1 50 100 150 200 250 300
TIME(SECS)
Figure 1.11: Output for low frequency sinusoidWiener model.
24
OUTPUT UNITS
.1
Figure 1.12': Output for low frequency sinusoidHammerstein model.
i
I
I
25
I
I
similar, the levels are different and the nonlinear and linear subsystems do not
commute implying that the two models are not equivalent. An analysis of this
I
i i
phenomenon is instructive, but it does not resolve the choice between the two
models and so is presented in Appendix A.
The model we will use for both subsystems of the muscle system is the Wiener
i
I
model. In fact, we can get either model to match human subject data by adjusting
the inputoutput curves associated with the nonlinearity. However, we feel that
the linear system forced by the work load followed by the nonlinear capacity
limiting function represents the physiologic reality better than the reverse order
of the Hammerstein model and so we will choose the Wiener model.
I
The linear subsystems will all be simulated as discussed in The Bicycle section.
However, here we do not assume that the input remains constant during the
internal time step and we approximate the integral in the evolution equation
(1.3) with the trapezoidal rule to get:
y(tj+1) = y{ij)exp{(At/Ti)) + (l/ri)(At/2)(H(tJ)ep((At/rl)) + U{tj+1))
! (1.4)
>1
where y is the output of the linear subsystem, At = tj+i tj, tj+i] is the jth
ll
internal time step and r, is the time constant for the respective subsystem. (In
the following Tqin will be the time constant for the muscle flow subsystem and
i
Tvm will be the time constant for the muscle oxygen consumption subsystem.)
'i
There are some clinical data regarding the actual values of Tqm and Tvrn ob
26
tained from isolated preparations of animal muscles and indirect observations in
man. There is also no a priori reason to assume that the off kinetics (response
to decreasing work load) are the same as the on kinetics. (See [23].) Barstow
and Mole [5] construct a model to represent VO2 under the influence of step
inputs. They use values of rqm from one to forty seconds and values of rvm from
ten to forty seconds and we will also consider values in this range. We will also
allow the on and the off time constants to be different. It is also possible
that there are actual delays in the system. (See [23].) These are modeled in the
simulation by simply shifting the time at which the input is applied. Of course,
all of these constants vary from one experimental subject to anotherjust as
j
each experimental subject has their unique set of parameters, so the simulation
corresponding to that subject will also have its unique set of parameters.
The next question to be answered is how to model the zeromemory, nonlinear
I
compartments. ,Note that when the input is constant the linear compartment
S
output approaches a steady state value as does the output of the cascaded model.
There are data available relating this steady state input to the steady state output
!
[31] and we will use these results to construct curves of both oxygen consumption
and muscle flowl verses input load for a general subject. We will then use these
curves as the nonlinear part of each subsystem. Again, these curves are unique
to each individual and must be determined for each experimental subject we
will do in Chapter Two. However, Figure 1.13 is a plot of the relation between
j
muscle oxygen consumption and apparent load presented by the output of the
27
Work load (watts)
Figure 1.13: Nonlinear relation between work load and muscle oxygen consump
tion.
linear dynamic subsystem constructed for a hypothetical subject from data in
several sources such as [31] and [10]. Figure 1.14 is a similar plot of muscle
blood flow verse;s apparent load. In each there are portions that are nearly linear,
especially if we'limit our considerations to small changes in the work load, and
each approaches a horizontal asymptote with increasing work load, as expected.
Note that we do not know the exact shape of these curvesthey are constructed
by interpolating and extrapolating data. Their shape suggests that they may fit
28
Flow (ml/sec)
I
Figure 1.14: Nonlinear relation between work load and muscle blood flow.
29
a MichaelisMenton curve:
i , ki input
, output = k3 + 7.
(2 + input)
There is some physiologic support for this in that we can view either subsystem as
flow through a limiting compartmentspatial limiting in the flow case and limits
on the rate at which oxygen can be consumed in the oxygen case. Similarly, the
MichaelisMenton relation attempts to model reaction kinetics as flow through a
system with limited enzyme available. At any rate, we are able to fit the empirical
I
data with the MichaelisMenton model and we will use them in the simulation
model.
It can be argued that the true purpose of the identification algorithm is to
determine the shape of these flow verses load and oxygen consumption verses
load curves for all load values. Suppose, for example, that we want to determine
V02max from submaximal inputs. Suppose further that the rate limiting factor
is flow. Then the input must be designed to determine, at least implicitly, how
flow extrapolates to its maximum. If the true shape of flow verses load in the
nonlinear subsystem were to be as in Figure 1.15, even the best designed test at
loads below 80 watts would give little hope of determining VO^max
Thus, in summary, we have modeled each muscle level subsystem as a Wiener
model. Figure 1.16 shows the result of a piecewise constant input into such a
model. The top plot shows the input of work load verses time. The second plot
shows how the linear system responds to this and finally the third plot shows
30
ANOTHER POSSIBLE INPUTOUTPUT RELATION
Figure 1.15: Another possible nonlinear relation between work load and muscle
blood flow.
31
INPUT TO LINEAR SYSTEM
Figure 1.16: Overall summary of Wiener model.
32
how the nonlinear system attenuates the output of the linear system to produce
the overall system output. This last plate shows both the output of the linear
system and the'final system output superimposed.
ii
There is one last step to complete the modeling of the muscle system. Recall
from Figure 1.5 that the Muscle System also has input from the arterial tree
I
and produces output into the venous system. The state variables that we will
observe in both [the arterial and venous sides is oxygen concentration in milliliters
of oxygen per milliliter of blood. Since we really want milliliters of oxygen to
1
be a measure of oxygen quantity, we could use a mass measurement or even
count molecules, but the tradition is to correct the volume of oxygen per unit of
blood to standard temperature and pressure and call the units milliliters percent.
I
From the above discussion we regard the local muscle blood flow and the local
muscle oxygen (consumption as a function of exercise load level. From this we
can determine the relation between oxygen concentration on the arterial inflow
I
side and oxygen concentration on the venous outflow side by a statement of
l
conservation of:!mass called Ficks principle. (See [40].) In this case it is:
; ((WO CWi)) Q{t) = Q02{t) (1.5)
where Q is thejlocal flow (ml/sec), Cao2 is the arterial concentration of oxygen
(ml oxygen/ml blood), C02 is the venous oxygen concentration (ml oxygen/ml
blood) and Q02 is the local oxygen consumption (ml oxygen/sec). It is known
that an exercising subject tries to maintain several arterial values in a narrow
33
range as long as the inspired gas mixture is not varied [40]. One such value
that is defended is arterial oxygen content at about 21 ml/100 ml blood [5]. We
have found by ear oximetry that this value does in fact decrease during vigorous
exercise. However, we will assume it to be a constant in the simulation and so
we can solve equation (1.5) for Cv02 as a function of local flow and local oxygen
consumption and, thus, of input work load. This then becomes an input to the
Venous System, which we model next.
The Vendus Subsystem
In view of the previous discussion, it is clear that as exercise load changes two
things will occur: Local flow into the venous system will change and the blood
that flows into the venous system will have varying oxygen content. We will view
the venous system as in Figure 1.17. There are two sources of inputone is flow
and oxygen content from the muscle box and the other is flow and oxygen content
from the visceral (or other compartment). (Recall Figure 1.5.) The output of
the venous system is the flow and blood oxygen concentration into the right side
of the heart. It is well known that the venous system can act as a variable size
storage reservoir and we will model this as three fixed volume pools and one
I
variable volume mixing chamber. Two of the fixed size pools connect the source
at the muscle box and the source at the visceral box with the mixing chamber
and the third fixed size pool connects the mixing chamber with the sink at the
right heart. We will we not take into account possible beat to beat blood storage
34
Figure 1.17: The model for the venous subsystem.
35
in the pulmonary vasculature and will require this flow into the right heart to
il
equal the left heart cardiac output.
The size of the venous pool has been determined by several methods and we
will use standard values for the equilibrium values and apportion them among
1'
the four pools in a manner that gives good agreement with experiment. The
']
input (flow and oxygen concentration) from the muscle compartment is a known
function of work load as derived in the previous section. The flow and oxygen
concentration from the visceral compartment will be determined by the methods
presented in the subsection titled The Visceral Compartment. Some previous
I
models of the cardiorespiratory system have assumed that blood from the muscle
i
compartment is! fully and instantaneously mixed with the venous pool. However,
we will track small elements as they exit the muscle compartment and traverse
the venous pool to the mixing chamber, treating this part of the venous system as
a first infirst out buffer. This part of the system between muscle and the mixing
chamber will be considered a constant volume pipe and each discrete bolus of
blood will have associated with it a content of oxygen that remains constant as
the element traverses this part of the venous system. The pipe will empty into the
I
i'
mixing chamber at a rate determined by the flow into the pipe, which is simply
]<
the output of the muscle compartment. The input from the visceral compartment
I
will be treated in a similar manner.
' t1
In the mixing chamber we will assume complete mixing occurs. Since we know
!
the initial size of the chamber and since we know the flow into it and flow out
:: 36
II
of it (cardiac output) we can determine the size of the mixing chamber at any

time and it will expand and contract depending on the various time constants
: 'l
i
of the flow systemscardiac output, muscle flow and visceral flow. Also we will
assume that the initial oxygen concentration in the chamber is the average rest
oxygen concentration of venous blood, .16 ml oxygen per ml blood. (See [31].)
Thus, since we Know the concentration of oxygen flowing into the chamber we can
'l
calculate the total oxygen content of the chamber and the oxygen concentration.
This oxygen concentration is then associated with each element of blood in the
mixing chamber at the current time. These elements then become the output
l
of the mixing chamber and we view them as being pumped out of the chamber
I;
each time step in a quantity equal to the cardiac output. We then track these
elements as they traverse the fixed volume pipe leading from the mixing chamber
l
to the right side of the heart and then into the lung and then into the left side
of the heart for ^pumping into the arterial tree.
I
Note that the transit time is flowdependent in that the time delay (TD) in
each part of the system is determined by an integral of the flow:
f flow dt volume
JtTD J
In the simulation we model this by tracking elements of blood through the venous
system as time evolves discretely. A continuous form of this model for the venous
system appears in the literature as early as the 1967 reference by Grodens [20]
and in a more recent model by Essfeld and Hoffman [14]. We feel it is much more
37
realistic than the complete mixing model. It also introduces a variable delay
nonlinearity into the system that is fascinating in its own right.
The Lung Subsystem
As we discussed briefly in the previous section we do not allow variable storage
of blood in the pulmonary vasculature. During each time step of the simulation
a quantity of blood enters and the same quantity of blood exits the lungs and
both of these equal the cardiac output. The input to this subsystem comes from
the fixed pipe connecting the venous mixing chamber with the right heart as
shown in Figure 1.17 and the output of this subsystem is the oxygenated blood
that is input into the left heart. This subsystem calculates the simulation output
variable that we are most interested in VOi Again we use Ficks principle and
make the assumption that the arterial oxygen concentration has a fixed value of
21 ml oxygen per 100 ml of blood. In its appropriate form Ficks principle is,
using the previous notation:
VO2 = (Ca02 ~ Cv02) Q
or, since arterial oxygen concentration = .21 ml oxygen/ml blood,
V02 = (21 Cv02) Q
Since we know the venous oxygen concentration and cardiac output we can cal
culate VO2 and we finally have this value as a function of work load input.
38
It is now possible to determine ventilation or other easily observable respira
i j
tory parameters by making assumptions about lung mechanics. We will not do
i
this since we are interested in VO2. However, we must note that although we
have oxygen consumption measurements available at time samples limited only
by the time discretization in the simulation, in reality one is practically limited
to observing the process only on a breathtobreath basis and the observation is
done at the mouth not at the alveoli. We will assume full mixing in the lungs
and assume a constant dead space. This makes our alveolar values the same
as the values collected at the mouth. We resolve the breathtobreath question
either by integrating our nearly instantaneous observations over a breath cycle to
produce breath to breath output or by assuming there is an interpolation algo
rithm, as discussed earlier, that essentially turns breathtobreath sampled data
into continuous data.
The Cardiac Subsystem
Cardiac output will be modeled in a manner similar to muscle blood flow.
We again use the Wiener model of a linear subsystem with first order dynamics
!
cascaded into a zero memory nonlinear subsystem. The input to the system is
work load and the output is the cardiac output, which as discussed before, is the
amount of blood pumped from the venous system through the lung system and
into the arterial system on each internal time interval of the simulation. More
complex models for the linear subsystem have been suggested (see Fugihara [16])
39
400
Work load (watts)
Figure 1.18: Nonlinear function between work load and cardiac output.
and since the simulation is modular they would be easy to implement. In this
thesis we will use only the first order system with time constant rco in the one to
forty second range and no fixed delay. Like the time constants for the subsystems
discussed earlier the values of rco vary from subject to subject. The nonlinear
l
subsystem is similar to the nonlinear subsystem for muscle flow. The known
steady state values from [31] or [10] are used to generate a table of cardiac outputs
as a function of, work load. As before, the intermediate values are determined
either by linear interpolation or by fitting to a MichaelisMenton equation. (See
Figure 1.18.)
40
The Visceral Subsystem
This part of the system is represented as 0 (for Other) in Figure 1.5. It
represents such body parts as viscera, brain, nonexercising muscle, skin, and
other physiologic compartments not otherwise represented in the model. We
model it by saying that the flow into it is the cardiac output minus the flow
i
to exercising muscle. We assume that its consumption of oxygen is relatively
fixed and independent of work load. At rest this compartment uses about 300
ml oxygen per minute or 5 ml oxygen per second [31]. The oxygen content at
the input to this compartment is the arterial value of .21 ml oxygen per ml
blood by the same argument used at the input to the muscle compartment. We
can, thus, apply Ficks principle again and determine the oxygen concentration
at the output. We assume no storage in this system and so the flow at the
output is the same as the flow at the input. This output flow and output oxygen
concentration become the input to the fixed venous compartment connecting the
visceral compartment to the mixing chamber.
The Arterial Tree
The final part of the system is the arterial tree that links the cardiac subsystem
with the visceral compartment and the muscle compartment in Figure 1.5. The
model we use is quite simple in that we will view it as a fixed volume pipe filled
with blood at the constant oxygen concentration of 21 ml oxygen per 100 ml
41
blood, as discussed earlier. This approximation allows us to essentially ignore
time delays and. flow on the arterial side.
The Overall System
Now that we have constructed each subblock of the overall system let us
again consider the overall system as an inputoutput system. As we have seen,
one possibility is to view the input as an element of Rn where the jth component
of this vector is the value of the input during the jth time interval. Each element
is held for one switch time to become the input to the subsystem downstream.
The inputoutput mapping in this case is from Rn to the output space.
Since we know the output of the computer (input to the bicycle) at all times
it is possible to, conceptually resample it at any rate that we choose including
one that is not equal to the time separating the components of the original input
vector. This is the view most commonly found in the physiology literature where,
for example, a total time of experiment might be 10 minutes (600 seconds). The
switching time could be taken to be 5 seconds (ie. the computer would switch
levels every 5 seconds) and so n = 600/5 = 120. (This is the n in the previous
paragraph.) The output of the sample and hold is then resampled at 1 Hz., thus
yielding one sample every second and so m = 600/1 = 600. We can, in this case,
view the experiment as a mapping from Rm to the output space.
A third possibility is to view the original R input as a sample from a finer
grained, even continuous, preinput. In this case the computer can be consid
42
I
ered a sampler in series with a zero order hold [15] and, in this case, the model
mapping is from a subspace of a function space limited only by the restrictions
imposed on its frequency content by the sampling theorem.
i
The output can be viewed in a similar way. One can accept the breathto
%
breath values actually output from the physiologic system, in which case the
I
output space is Rfc, where k is the number of samples taken of the output. Or
one can smooth] them as discussed in the section about the measuring devices.
If this smoothed output is used directly, the output space is a function space. If
the smoothed output is resampled, the output space is then Rp, where p is the
number of resamples. These possibilities will become important in Chapter Three
where we will apply methods to identify systems from inputoutput data and they
will be discussed further there. Also, since it is really the subject subsystem in
Figure 1.3 that we are interested in, and since we know the transfer functions for
I
the bicycle subsystem and the measuring devices subsystem we can invert these
{
subsystems and [study just the mapping from the input to the subject subsystem
to its output. Along these same lines we will consider the question of what input
signal can be input to the overall system to produce a desired testing signal at
;i
the input to thef subject subsystem.
[ An Example
I
In this section we will present an example of how the simulation works. Since
we feel that the 'model represents the physiology we will not add the word simu
43
'i
i
Figure 1.19: Work load in watts that is input to the example system.
lated in the following and, for example, say muscle blood flow for both muscle
blood flow in the human experimental system and the simulated muscle blood
!
flow in the simulated system. Figure 1.19 shows a sample input choosen at ran
dom, subject only to the requirement that its amplitude range is limited to the
physiologic range of zero watts to about 250 watts. In the real world experiment,
'I
the investigator decides how frequently to apply the sample and hold to the signal
(say every ten seconds as we will use here) and he enters these corresponding val
j
ues into the computer. When the experiment is started these values are clocked
out to the bicycle through the D/A converter. Figure 1.20 shows this output. As
;l
we discussed earlier we will process the data digitally and so we must pass the
sampled and held data through an antialiasing filter. In this case we will sample
at 1 Hz. and so the Nyquist frequency is 0.5 Hz. Figure 1.21 shows the result
after the antialiasing filter. Since few high frequencies were present in the original
signal and we are sampling at a relatively high frequency we see little effect from
44
Figure 1.20: Original input (solid) and the sampled and held signal (dashed).
Figure 1.21: Output of the sample and hold after it is passed through the an
tialiasing filter. '
the filter except! a hint of the Gibbs phenomenon. We assume that the bicycle
subsystem is linear with time constant = 1 sec. The output of the bike, which
is the load the subject actually rides against is shown in Figure 1.22.
This input acts on the three Wiener subsystems and produces the correspond
I
ing outputs. Figure 1.23 shows the resulting muscle blood flow verses time. If we
compare this output with the original input we see a time lag as we would expect
and we can also; see the saturation effect of the nonlinearity. Figure 1.24 shows
45
UNITS UNITS
200
TIME(SECS)
Figure 1.22: Output of the bicycle subsystem.
Figure 1.23: Muscle blood flow in ml. per second resulting from the input shown
in Figure 1.19.
i
46
Figure 1.24: Muscle oxygen consumption in ml. oxygen per second resulting from
the input of Figure 1.19.
the resulting muscle oxygen consumption verses time. Note that the muscle blood
flow and the muscle oxygen consumption curves are similar. This follows since
the time constants for the two systems are chosen to be the same in this example
and since the curves representing the nonlinearities are also similar. Figure 1.25
(top plate) shows the resulting cardiac output verses time plot from the input in
Figure 1.19. Again it is similar to the two previous plots. The bottom part of
Figure 1.25 shows the output of the linear subsystem of the Wiener subsystem
for cardiac output. Comparing the overall output (top) with the intermediate
output of the linear system we see the effect of the modulating nonlinearity.
The muscle blood flow and muscle oxygen consumption outputs combine as
discussed earlier to produce the oxygen concentration of blood entering the venous
pool and a plot of this quantity is shown in the top part of Figure 1.26. The
corresponding concentration of oxygen in the venous pool is shown in the bottom
47
I
Figure 1.25: Cardiac output resulting from the input in Figure 1.19 (top) and
the corresponding output from the linear subsystem of the cardiac output Wiener
system (bottom).
!
i
li
48
Figure 1.26: Oxygen concentration of blood entering the venous pool (top) and
the oxygen concentration of blood in the venous pool (bottom) resulting from
the input in Figure 1.19. Both graphs have units of ml. oxygen per ml. blood.
49
Figure 1.27: Oxygen consumption of the system in ml. oxygen per second result
ing from the input work load in Figure 1.19.
plot of Figure 1.26. Note that the result is smoother than the top plot and lags
behind it, both j a result of the fact that the mixing chamber is an integrating
device. The blood is then pumped from the mixing chamber and past the lungs
where it takes up oxygen. This determines the final oxygen consumption and,
thus, the final output of the system and is plotted in Figure 1.27.
Conclusion
At this point we have developed a computer simulation model of the human
.1
cardiorespiratory system and an associated experimental environment that takes
as input a collection of numbers that represent work load and an associated
collection of times to switch from one work load to another. The model produces
an output consisting of observations of VO2 and several internal variablessuch
as cardiac output, venous oxygen concentration, and venous volume. If we assume
I
the dynamics sire first order we have fifteen parameters. Twelve of these are
represented as an on and an off time delay and an on and an off time
constant for each of the three linear subsystems discussed above. We use on
here to mean that the input to the system is increasing and off to mean that
tl
h
the input to the system is decreasing. The three other parameters are the sizes
of the fixed venous compartments.
We also have initial conditions for each of the systems and the initial size of
i
the mixing chamber. The three curves of steady state values verses work load
are also variables in the simulation. The system has nonlinearities that arise
from several sources. One is the application of Ficks principle; for example,
the calculation that determines venous oxygen content, requires the quotient of
oxygen consumption and blood flow and so introduces a quotient nonlinearity.
Another application of Ficks principle is in the calculation of oxygen consumption
at the lung. Here it requires the product of venous oxygen content and flow and so
introduces a product nonlinearity. Other nonlinearities arise from the nonlinear
portion of the Wiener model used to model cardiac output, muscle blood flow and
i
muscle oxygen consumption. A third nonlinearity arises from the complicated,
time delay that occurs in the venous system.
We believe the idea of using a Wiener model with the nonlinearity represented
as we have done here is new to cardiorespiratory modeling. It allows for com
i
plicated inputs and also allows different on and off dynamics. It also establishes
that an important part of the identification of the system is the determination of
' 51
the inputoutput curves for the nonlinear subsystems. We also believe that the
venous model where each bolus of blood is tracked across the venous pool is also
new.
In the next chapter we will present and discuss sample experimental runs and
apply some techniques from systems analysis to the simulation. In Chapter Four
we will use the simulation to apply and compare time domain system identifica
tion techniques and show how these results have application in the actual testing
1
of the human subject.
t
52
CHAPTER 2
ANALYSIS OF THE MODEL
Introduction
In this chapter we will determine the output of the simulation under several
inputs often used in system analysis ramps, steps, sine waves and a class of
signals often called constant switching pace (CSP) [30] that are useful in testing
and identifying systems. We will discuss why the simulation responds as it does
and compare the simulated output with output from human subjects. In the last
section of this chapter we will introduce a noise model into the simulation.
!
I
Ramp Inputs
it
The first input that we will consider is a ramp as shown in Figure 2.1. This
# l] *
is the type of input used clinically to determine Vr02max In such an experiment
the braking voltage to the bike is slowly increased (and so the work load) until
the VO2 plateaus, the subject fatigues or stops for some other reason. Figure 2.2
shows an actual1 response from a human subject. In this plot we have averaged
the breathbybreath output over eight consecutive breaths using a standard pro
tocol used in exercise testing [24] and so this output appears relatively smooth.
I
Figure 2.3 shows the human response with the simulation response superim
OXYGEN UPTAKE(ML/MIN)
!
)
I
Figure 2.1: Sample ramp input of work load vs time.
Figure 2.2: Output of a human subject under the ramp input.
54
o
w
in
\
S
E
D
CL,
H
O
I,
Figure 2.3: Human subject output (dashed) and simulation output (solid) under
i
a ramp input. '
posed. They agree quite well and, in fact, we can make them agree very closely
i'
because we use this test to fit each of the cardiac outputwork load, muscle oxy
gen consumpticjinwork load and muscle blood flowwork load curves discussed
in Chapter One to the given human subject. These three curves determine the
V02 at each input level and we could determine them by testing the subject
with many different constant inputs. Instead, we choose the three curves above
by trial and error in such a manner that they produce output that matches the
ramp. Note that since there are three curves at each data point we have three
degrees of freedomone for each curve. This process only uses up one degree of
freedom at each point and so we have two left to fit further subtleties in the data
J
if needed. !
!'
One thing to note is that the output is nearly linear over a large input range.
This is sometimes cited as evidence that the V02 system is linear. In reality, this
il
55
apparent linearity is an artifact of the experiment. Recall from Ficks principle
that the oxygen consumption at the lung is determined by the cardiac output
jl
!,
and the oxygen! concentration gradient between the blood flowing into the lung
'i
and the blood flowing out of the lung (recall Chapter One):
V02 = (.21 CV02) CO,
where .21 is the oxygen concentration of the arterial blood, CV02 is the oxygen
il
concentration of the venous blood and CO is the cardiac output. Under the
,1
influence of a ramp input we see from Figure 1.18 that cardiac output rises and
then begins to saturate as load increases. We also know that as load increases ve
1
nous oxygen concentration falls as the exercising muscles extract more and more
oxygen from the blood flowing through them. Thus, as the cardiac output (CO
above) saturates, the venous oxygen concentration (Cvo2) falls. These comple
ment each othet in the above equation to produce the nearly linear output until
at a high input level, the effect breaks down and the limiting effects of both C0
and Cvo2 appear and VO2 plateaus.
This complementary effect does not necessarily hold with other inputs and
1
shows how important it is to use a variety of inputs in testing the system. In
fact, what we see here is that the input prepares the venous pool in such a
way that its decreasing oxygen concentration complements the increasing cardiac
output to produce a nearly linear system output. Consider instead inputs that
prepare the venous pool in different ways as shown in the top plate of Figure 2.4.
TIME(SECS)
i
Figure 2.4: Step! up and step down inputs (top plate) and the resulting simulation
outputs (bottom plates).
: i
The solid line is a step down and prepares the venous pool to be oxygen depleted.
The dashed line'is a step up and prepares the venous pool to be oxygen rich. The
resulting outputs are shown in the bottom plate of Figure 2.4. If the system were
truly linear the sum of the outputs would be the output of the sums (this is the
definition of linearity). The top plate of Figure 2.5 shows the output resulting
from the sum of; the inputs and the bottom plate shows the sum of the outputs.
>1
Clearly they arei different and the system is not linear.
There are several reasons why the human VO2 vs load system is not linear.
I
'I
57
I
Figure 2.5: The output of the simulation under the summed step inputs (top
plate) and the; sum of the outputs under each step input separately (bottom
plate).
i
I
!
58
100 
co
H
Ã‚Â£
50 100 150
TIME(SECS)
Figure 2.6: Sample step input.
I
One is the complicated dynamics of the venous pool as discussed above. Another
!
is that the time constants in the true human system are not exactly the same
for the off dynamics as for the on dynamics and a third is that we allow the
possibility of pushing the subject into the plateau portion of the curve in Figure
2.2 and this clearly produces nonlinear results. Further experimental results
showing that the human system is nonlinear and a discussion can be found in
[23].
Step Inputs
The next input form we consider is the step input. A sample input is given
in Figure 2.6 with the corresponding output in Figure 2.7. By averaging the
output over many repetitions of real experimental data and by fitting curves to
the data several investigators (for example, [23]) find that the human subject
has the biphasic output of Figure 2.7. (Note in Figure 2.7 that for about fifteen
59
1 Figure 2.7: Simulation under the step input.
seconds after the step in the input that the shape of the output is different than
the later portion giving the output a biphasic appearence.) Whipp et al. [46]
speculate that the early phase consists of a rapid increase in cardiac output which
is emptying a venous pool prepared by the low load portion of the step input to
i
be oxygen rich. The second phase occurs as the venous pool becomes more and
more oxygen depleted as a result of receiving oxygen depleted blood form the
muscles responding to the high load portion of the step. Under a step input our
simulation becomes equivalent to the model proposed by Barstow and Mole [5].
They show that their model is in agreement with experimental step data and so
it is not surprising that our model also matches step data well. Figure 2.8 shows
the simulation output with the output of an experimental subject superimposed.
They agree well. In fact, we use this step data to determine the time constants of
the linear portions of the underlying Wiener subsystems of Chapter One. Again,
there are three degrees of freedomone for the linear subsystem in each of the
60
30
Figure 2.8: Output of the simulation (solid) with output of a human subject
(dashed) superimposed. Output units are ml. oxygen per second.
!
three Wiener systems. We find that we can fit the data relatively well by choosing
the three time constants to be the same and so we only use one degree of freedom
in this fitting. This leaves two to model further subtleties in the data.
To see how the choice of these time constants can effect the fine structure
of the model consider the relation between the time constant for the Wiener
subsystem that determines muscle blood flow (r?m) and the time constant for
the Wiener subsystem that determines muscle oxygen consumption (rm). If
Tgm < Tvm an increase in load causes muscle blood flow to increase slower than
muscle oxygen consumption and the effluent from the muscle falls in oxygen
concentration as one would expect. However, if rgm > rm, flow would increase
faster than oxygen consumption and we would have the paradoxical result that the
oxygen concentration of blood flowing out of the muscle would actually rise with
increasing load.: By varying the relative magnitudes of these two time constants
61
150
co
U 100
50
_______i_______i______i_____i________i_____i________i_____i________i______
0 20 40 60 80 100 120 140 160 180 200
TIME(SECS)
Figure 2.9: Soule responses that can be obtained by fixing Tqm and varying rvm.
we can adjust the fine structure of the response to fit the data for a given subject.
Figure 2.9 shows some of the variety of responses we can produce.
Sine Wave Input
The next input we consider is a sine wave. Figure 2.10 shows a .01 Hz. sine
wave in the time domain in the top plate and in the frequency domain in the
bottom plate. When this signal is sampled and held we get the more complicated
signal shown in the time domain in Figure 2.11. In this example we will sample
at ten second intervals. It is known [25] that sampling in the time domain has the
62
I
1 :ri11ir
O !
5 
!
o '
CL,
o .______________,__________i_______i_______,_________,,________
0.2 t0.15 0.1 0.05 0 0.05 0.1 0.15 0.2
, FREQUENCY(HZ)
Figure 2.10: Sine wave in the time domain (top) and frequency domain (bottom).
Frequency spikes are at 0.01 Hz. and 0.01 Hz.
(
Figure 2.11: Sampled and held sinusoid in the time domain.
63
PSD
Figure 2.12: Spectrum of the sinusoid after sampling,
effect of replication of the spectrum in the frequency domain shifted successively
by A
sinusoid will have the frequency domain representation shown in Figure 2.12.
This sampled signal is then held constant for ten seconds until the next sample.
This function is called a zero order hold (ZOH) [15] as opposed to a first order
.1
hold where two! successive samples are linearly interpolated or a second order hold
where successive sample are fit by a quadratic or other higher order holds. In the
i
time domain this function can be represented as in Figure 2.13 which shows the
impulse response of the ZOH. The transfer function of this system is the Fourier
transform of the impulse response and so is just the Fourier transform of a shifted 1
1This argument is often presented in conjunction with analyses of aliasing, where it is pointed
out that if the samples in the time domain are too far apart (T too large) then Au will be too
small and the shifted copies of the spectrum will overlap resulting in aliasing.
64
INPUT TO ZOH
, (IMPULSE)
OUTPUT OF ZOH
(STEP)
ZOH
TIME
TIME
Figure 2.13: Effect of ZOH on an impulse at t = 0.
rectangular window:
\
! Hzoh = exp(jwT/2)Tsinc(u>T/2),
where ainc(x) = sm(x) f if ^ 0 and sinc(x) = 1 if x = 0. In terms of the power
transfer function, we have:
H*zoh{3u)Hzoh{jw) = T2ainc2(uT/ 2),
This is plotted in Figure 2.14. We can now cascade the sampling action with the
I
hold action on the sinusoid by multiplication of the power transfer functions in
the frequency domain to determine the output of the ZOH. Figure 2.15 shows
this calculation jby superimposing the graphs. From the point of view of testing
i
the system the important point is that we are not presenting a sinusoid to The
Bicycle subsystem and, thus, to The Subject subsystem, but a more complicated
signal consisting of several frequency components. We would prefer to input a
65
POWER UNITS
FREQUENCY
Figure 2.14: Power transfer function for the ZOH in normalized frequency units
of (wT/2)..
66
POWER UNITS
FREQUENCY
Figure 2.15: PSD of the sampled sinusoid superimposed on the power transfer
function of the ZOH. Frequency is in units of u>T/2.
67
cutoff frequency of 0.1 Hz.
sinusoid into the human system and the above suggests that we might do this by
filtering. Recall that we still have the antialiasing filter to impose before passing
the signal to the bike. If we were using this filter only to prevent aliasing the
cutoff would be near the Nyquist frequency. If we are resampling (recall Chapter
One) at, say, one sample every 0.5 seconds, the Nyquist frequency would be
//v = 1Hz. However, we can certainly set it lower. Figure 2.16 is a plot of the
output of the antialiasing filter if the cutoff is set at 0.1 Hz. If we set the cutoff
frequency even lower we can finally remove all the frequency components except
the fundamental as shown in Figure 2.17. Figure 2.18 shows both the original
I
signal input to the computer and the resulting filtered output that is then input
to the bike. We see that they are very similar except for a phase shift. Also
note that the bicycle itself is a low pass filter. Recall from Chapter One that we
model it as a linear system. Figure 2.19 shows the power transfer function for a
68
l
Figure 2.17: Sampled and held sinusoid after filtering through low pass filter with
cutoff frequency of 0.05 Hz.
Figure 2.18: Input to the computer (solid) and final output to the bicycle
(dashed). 
t
I
69
1,
Figure 2.19: Power transfer function for a typical bicycle.
typical bicycle and Figure 2.20 shows the result of passing the sinusoid through
the bicycle alone without the antialiasing filter. It is clearly nearly sinusoidal.
i
Thus, by proper filtering we can test the human system with a sinusoid in spite
I;
of the obligate sampling and holding that are part of the system.
Figure 2.21 Jshows the result of testing a human with a 0.01 Hz. sinusoidal
input with the output of the simulation superimposed. Clearly the two agree
very well. Here we stress that in this case the simulation is not fitted to this
i,
j
data. The amplitude response of the simulation is determined by fitting the
70
I
50
, TIME(SECS)
i
Figure 2.21: Output of simulation and human subject under sinusoidal input.
j
ramp response and the time constants of the simulation are determined by fitting
the step response as discussed earlier. There is no fitting to the sine response
yet they agree very well. The mean error (the mean of the quantity determined
by subtracting the output predicted by the model from the observed value) for
this study is 0.52 units. After we remove this bias the standard deviation of the
residuals is 2.52 units. We will see later in this chapter that the noise has a
standard deviation of 2.6 units about its mean and so this represents very good
agreement. Figure 2.22 is the autocorrelation function of the residuals showing
that they are nearly white. Figure 2.23 is the crosscorrelation function of the
r
residuals with the input showing little crosscorrelation. All of these factors show
that the simulation is doing a very good job of modeling the subject for the
sinusoidal input data.
71
1
3
0.5 
0 
0.51''1'1111
0 2 4 6 B 10 12 14 16 18 20
LAG(SECS)
Figure 2.22:, Autocorrelation of the residuals for the sinusoidal input case.
1
0.5
OT
M o
0.5
1
0 5 10 15 20 25 30 35 40 45 50
LAG(SECS)
Figure 2.23: Crosscorrelation of the residuals and the input of the sinusoidal
input case.
72
'I
Constant Switching Pace Inputs
In this section we consider the output of the simulation and the human under
the input of a constant switching pace (CSP) signal. We will discuss this signal
class more fully later and use them in the identification technique that will be
i
developed. The signal is generated by selecting a hold time or viewed another
way the time between switches from one level to another. During each hold
time we choose the input level by a random draw from a uniform distribution.
I
\
Figure 2.24 shows a CSP input and the resulting output of both the simulation
l
and the human system. Again the outputs agree very well and as in the sinusoidal
input case the simulation is not fitted to this specific signal. In this experiment
the mean error was .43 and the standard deviation of the residuals was 2.69.
l The Model with Noise
In this section we will introduce a simple noise model into the simulation.
I,
Figure 2.25 shows a time domain plot of output from a human subject about its
I
mean in the top plate and a frequency domain plot in the bottom plate under a
constant input., The standard deviation of this signal is 2.6 units. Figure 2.26
shows a histogram of the amplitude distribution of this time series in the top plate
and its autocorrelation function in the bottom plate. The histogram appears
normally distributed. To verify this we plot in Figure 2.27 the residuals against
73
I
1 TIME(SECS)
I
Figure 2.24: Constant switching pace input (top plate) and the corresponding
output of the human system (dashed line, bottom plate) and the simulation (solid
line, bottom plate).
74
PSD UNITS UNITS
10,
100
200 300 400
TIME(SECS)
500
600
0.5 
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
FREQ(HZ)
Figure 2.25: Time domain plot of the output from a human subject about its
mean (top plate) and the corresponding frequency domain plot (bottom plate).
75
AUTOCOR. UNITS
BIN
Figure 2.26: Histogram of the noise amplitudes (top) and autocorrelation function
of the noise (bottom).
76
4
VALUE
Figure 2.27: Plot of residuals against the normal scores for noise data,
the normal scores calculated by the MINITAB 2 routine NSCORE [33]. The
plot is nearly linear as one would expect if the data were indeed normal. When
we calculate the correlation coefficient of these data with a straight line we get
.988 which allows us to accept the hypothesis that the data comes from a normal
distribution by the ShapiroWilks test [33]. From the autocorrelation function
of Figure 2.26 we see that the noise is nearly white. Thus, at least to a first
approximation the output noise is nearly white and nearly Gaussian. The near
Gaussianity can be explained by assuming that the final noise results as the
sum of many small, independent sources and then applying the central limit
theorem to see that the final noise is near normal. In order to understand the near
whiteness of the noise we observe that the cascaded subsystems are smoothing
and so noise effects tend to be removed as data flows through the system. Each
Wiener subsystem consists of a linear system which acts as a low pass filter.
Registered trademark of Minitab, Inc.
77
200
Figure 2.28: Noisefree input (top) and corresponding output (bottom).
Similarly, the venous pool acts as a smoothing system. It averages the inputs
with the previous state of the pool to produce the current oxygen concentration
in the pool. All of these effects tend to remove the noise effects from previous
stages. This is illustrated in Figure 2.28 which shows the result (bottom plate)
of inputing a noisefree signal (top plate) into the simulation. Figure 2.29 shows
the same input (top plate) corrupted with Gaussian white input noise with power
level at ten percent of the original input. The output, shown in the bottom plate,
i
is essentially unchanged from the output in the noisefree case. Thus, we will in
this thesis use additive, white Gaussian output noise as our first approximation
78
OUTPUT(ML/SEC) OUTPUT(ML/SEC)
300
Figure 2.29: Noisy input (top) and corresponding output (bottom).
79
to a noise model in the simulation. As more data on the exact values of the noise
in each subsystem become available we will incorporate it.
Summary
In this chapter we have studied the output of the simulation and the output of
human subjects under several standard inputs and introduced a noise model into
the simulation. The results show very good agreement between the simulation and
experimental data. An exhaustive validation would require fitting the simulation
to many subjects under a variety of experimental conditions. We simply do
not have the resources available to do this and will leave this for future work.
However, we feel that the physiology we have used to develop the simulation
is correct and that the examples presented here show that the simulation is an
excellent first order approximation to the human system. We also feel that its
incorporation of simple and complex nonlinear behaviors makes it very useful as
a testbed for the development of system identification algorithms. In the next
chapter we will develop such an algorithm and in Chapter Four we will apply it
to the simulation.
80
CHAPTER 3
DEVELOPMENT OF THE IDENTIFICATION METHOD
Introduction
In this chapter we will begin the process of constructing a method of system
identification. We here again state our objective of determining, by using only
input and output data, a global identification in terms of polynomial functions
of lagged inputs and outputs. If we were dealing with human subjects the objec
tive would be to design an experimental protocol that would be run on a given
subject to produce a mathematical model of the subject. Such a mathematical
model would be ideal if it produced the same output as the subject for every
possible input. Practically, such an ideal model can not be constructed and we
instead construct models that approximately match the subject output for some
limited class of input signals. In this thesis we will instead be identifying the
computer model of the human subject developed in Chapter One. Although the
identification of this relatively complicated model is interesting in its own right,
our idea is to use it to determine and solve problems associated with identifying
the human system. The problems that arise are parallel in both systems; how
ever, the simulated system allows us much more control. For example, in the
computer simulation we know exactly the structure of the system, its parameters
81
and its noise sources. We also have access to internal variables not available in
human experiments and we have the ability to alter all of these quickly and easily
without putting human subjects at risk or exposing them to long, difficult exper
imental sessions. Because the techniques we develop relate to both the simulated
and real systems, we will freely mix notation in the following and say, for exam
ple, the input to the bicycle is passed to the subject when, in fact, we mean
the input to the simulated bicycle is passed to the simulated subject. If there
are situations where this association of human based experiment and simulation
based experiment break down we will point them out.
In developing the identification there are several questions that must be ad
dressed. According to Ljung [29] we must determine three entitiesthe data
record, the set of candidate models and the rule for assessing the models. In
cluded in the data record are such questions as the form of the input, what
outputs to measure, when to sample both the input and the output and other
questions related to data collection. The candidate models are determined by
how we represent the system and the assessment rule is the criterion used to de
termine which model in the candidate set best represents the system. We must
also establish an identification method that selects from the candidate models the
ones that optimize the assessment rule. In the first part of this chapter we will
discuss each of these aspects in turn. In the second part of the chapter we will
apply them to the identification of a simple system in preparation for applying
them to the full simulated system in Chapter Four.
82
The Data Record
The data record we consider is an inputoutput pair where the input is the
test input the investigator enters into the computer and the output is an interpo
lated and smoothed version of the breathtobreath oxygen consumption data as
discussed in Chapter One. Regarding the input data, there are two constraints
imposed by the experimental system. One is that the input level is limited so
that the work load presented to the subject is between zero watts and about
250 watts. While it would be possible to design a driven bicycle and a coupling
mechanism so that energy transfer would be from bike to subject, this is clearly
a different system than we are trying to test and so the lower limit is unloaded
'I
cyclingzero watts. At the other extreme, there is a limit above which the subject
cannot maintain constant pedaling speed on the bike and so we must limit the
upper range of the input to about 250 watts. The second constraint is that the
l
duration of the experiment must be no longer than about thirty minutes. This is
because subjects fatigue and if the experiment is too long physiologic parameters
change and the system is no longer stationary. Both of the above constraints are
subject dependentan elderly patient with emphysema will have different limits
than a young endurance athlete. The ranges considered here are for healthy sub
jects and we will discuss the effects of limiting both load range and experiment
duration.
Another limitation on the input arises from the experimental setup as dis
83
Figure 3.1: The general modeling scheme.
cussed in Chapter One. Although it would be possible to drive the bicycle with a
continuous analog signal, the typical experimental setup consists of sampling the
continuous input at a constant rate and presenting this sample to the bike for a
constant hold time equal to the interval between samples. This technique of sam
pling followed by a zero order hold is an artifact of the experimental hardware,
but is also a common technique in control theory [2]. Figure 3.1 shows the over
ii
all modeling scheme. Arbitrary continuous inputs (A) pass through the Ideal
system to produce the corresponding continuous outputs (a). In reality, since
the experimental system can pass only piecewise constant inputs to the bicycle,
the actual experimental loop is (A) to (B) to (b) to (a). Thus, the arbitrary in
84
puts at (A) pass through the Sample and Hold to become the piecewise constant
inputs at (B). These then are presented to the True System. Inside the True
System these inputs produce the breathtobreath data that are interpolated
as discussed in Chapter One to produce the continuous outputs at (b). These
outputs are equivalent to the output of the Ideal System in the sense that the
distance between them in the chosen norm is less than some specified tolerance.
It is this tolerance that will determine the upper bound on the hold time: If the
hold time is too long, the signals at (B) will be too different from the signals at
(A) and the resulting signals at (a) and (b) will differ by more than the required
tolerance.
Figure 3.2 shows a .01 Hz. sine wave and the resulting sample and hold signal
for hold times of 5 seconds (top plate) and 20 seconds (bottom plate). It is clear
that the 20 second hold signal is less like the original signal than the 5 second
hold version. Therefore, we would expect that a 5 second hold in the sample
and hold of Figure 3.1 would produce a signal at point (a) more like the result
produced by the Ideal System than would a 20 second hold. The question of
whether they are within the specified tolerance is one we will quantify later in
this chapter when we consider a specific system.
There are many different possibilities for determining the hold time structure
in the identification experiment. One method is to let it be a random variable
as discussed in Tulleken [43]. In this thesis we will require that the hold time be
a constant. The determination of this constant is based on the general tradeoff
85
UNITS UNITS
Figure 3.2: A .01 Hz. sine wave and its sample and hold for hold times of 5
seconds (upper) and 20 seconds (lower).
86
that if the hold time is too long, then signals at points (a) and (b) in Figure 3.1
do not agree within the tolerance discussed above. On the other hand, if the
hold time is too short then the input will be deficient in low frequencies and
will not stimulate important modes of the system and will not allow adequate
excursion in the amplitude of the signals input to the nonlinearities to allow
complete identification.
The other factor we must consider is the amplitude during each hold time. If
the input is a sine wave or other deterministic signal the amplitude during the
hold time is the value of the input at the sample time. In designing test inputs
we will take a slightly different point of view (recall the discussion in The Overall
System section of Chapter One) and choose the amplitude value during a given
hold time as a draw from a probability distribution. Several possibilities have
been studied. Marmarelis [30] considers, among other possibilities, independent,
Gaussian values. Barker [4] and many others consider values that alternate be
tween two values according to a very specific algorithm which produces binary
pseudorandom values. In this thesis we will follow Billings and Voon [8] and use
independent, uniformly distributed values. In the analysis of a linear system or
in the local analysis of nonlinear systems there is much to be said for pseudo
random sequences. By construction they have peaked autocorrelation functions
(and therefore broad power spectral densities) that are very useful in identifica
tion [13]. Their shortcoming in the global identification of nonlinear systems is
illustrated in Figure 3.3. Each curve represents the inputoutput function of a
87
Figure 3.3: Two nonlinear inputoutput functions.
zeromemory nonlinear system. This means that for a given input on the xaxis
the output is the corresponding y value. The two curves have the same y values
at the x values of 25 and 75. If we were to choose these as the values for the
two levels the systems would appear the same. Just as a line is determined by
two points such bilevel methods work for linear systems, however, for nonlinear
systems they do not work in general. Gaussian distributed values also have a
similar shortcoming as illustrated in Figure 3.4 which shows a zero memory non
linear transfer function with two Gaussian distributions superimposed. One way
to view the identification process is that we approximate the nonlinearity with
polynomial terms in such a way as to minimize errors at points along the xaxis
88
OUTPUT UNITS
Figure 3.4: Nonlinear inputoutput function with superimposed gaussian distri
butions.
89
drawn from the distributions. Clearly the coefficients of the polynomials and,
therefore, the model resulting from identification will be different for the two
different sampling distributions. Since we hope to find a model valid across the
entire range of amplitudes the uniform distribution is an obvious choice. In the
nonlinear systems we consider, which are cascades of linear and zeromemory
nonlinear systems, there are other ways to insure adequate excursions in ampli
tude at the input to the non linearity, but in this thesis we will use the uniform
distribution and the amplitude values during the hold times will be uniformly
distributed across the range of the input. Thus, the test inputs that we will use
for creating the data record will be signals that switch values at multiples of a
fundamental time called the hold time. During each hold time the value will be
determined by a sample from a uniform distribution between the lowest possible
input value and the highest possible input value. An example of such a signal is
given in the top plate of Figure 3.5 for a 5 second hold time. Such signals are
often called staircase signals because of their appearance, but we will call them
constant switching pace (CSP) signals [30].
In the actual identification we will use methods that involve the cross and
autocovariance properties of the inputs and outputs. One useful property is
for the autocovariance function of the input to approach a delta function. A
generalization of this is a concept which deals with the condition number and
invertability of the matrix of lagged autocorrelation functions called persistent
excitation (see Ljung [29]).
90
Figure 3.5: Constant switching pace signal with 5 second hold time (top plate)
and its autocovariance function (bottom plate).
It is easy to show [41] that the autocovariance function for the inputs we are
using (CSP signals) is
r(r) =
where N is the hold time and that our input is persistently exciting of all orders.
This autocovariance function is plotted in Figure 3.5 (bottom plate) for N=5 and
consists of a triangular pulse. It is well known that as such functions become
0
N
0, otherwise.
more localized in time their Fourier transforms spread in frequency. This will
91
allow us to match the bandwidth of the input to the bandwidth of the system as
we will discuss in detail when we consider the specific system.
There are two further points we must consider. One is that the experimental
conditions require that the identification be done in a limited time and so we
are not really dealing with infinite sequences. This raises the question of how
many samples are needed for the finite sequences to approximate the properties
of the infinite sequences closely enough to be of use and the related question
of determining if some finite realizations are better than others in the practical
identification. The other point is that in linear system identification we need only
know the properties of the inputs and outputs up to their second momentsie.
means and covariances or corresponding spectra. In nonlinear systems we must,
in general, know higher order moments or polyspectra for a complete analysis
[35]. This raises the question of the higher order properties of the inputs we are
using. Neither of the above points appear directly in the identification method we
will develop and so we will not discuss them further in this thesis. The question
of the best finite realization, however, is quite important practically and we will
include it as future work.
To complete the discussion of the data record consider the sampling of the
output during the identification experiment. It this thesis we will sample the
output synchronously with the input. There are many other possibilities as dis
cussed in Crochiere and Rabiner [12]. In nonlinear systems there may be some
advantage to sampling the output faster than the input because harmonics of
92
the input frequencies can appear. We choose, however, to keep the processing
bandwidth the same for both input and output and deal with this problem by
widening the overall processing bandwidth. This will be discussed further in the
setting of the specific systems.
%
Candidate Models
There are several possible candidate mathematical models we could use [1], [3].
There are also several ways to deal with the problem of identifying a naturally
continuous system [37], [44]. However, in this thesis we will sample both the
input and output processes and construct a discrete model between these sampled
processes. Once this discrete model is identified, the way the overall model will
work is by sampling the given continuous input, passing these samples through
the discrete model to produce output samples and then these will be interpolated
to produce the final continuous output. This procedure is shown by the lower
loop in Figure 3.1. Arbitrary continuous inputs (A) pass through the Ideal
System to produce the corresponding continuous outputs (a). The corresponding
model loop will be (A) to (B) to (C) to (c) to (b) to (a). Although this process
may appear complicated it is commonly used in both system identification and
control applications. The principal advantage being that the discrete model can
be determined using digital computing techniques.
In this thesis we will restrict ourselves to two general types of models. One
is the Volterra series as presented in Rugh [36] and the second is the Nonlinear
93

PAGE 1
A METHOD OF SYSTEM IDENTIFICATION APPLIED TO A SIMULATION OF HUMAN EXERCISE by Hugh Austin King B.S., Iowa State University, 1968 M.S., New York University, 1970 M.D., University of Pennsylvania, 1977 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Doctor of Philosophy Applied Mathematics 1992
PAGE 2
This thesis for the Doctor. _of Philosophy degree by Hugh Austin King has been approved for the Department of Mathematics by William Briggs David Fox Harvey Green erg Richard Jones 4/3D ;cp Date
PAGE 3
King, Hugh Austin (Ph.D., Applied Mathematics) A Method of System Identification Applied to a Simulation of Human Exercise Thesis directed by Associate Professor William Briggs A simulation model of the inputoutput system consisting of a human riding an exercise bicycle is developed. The input consists of work load and the output is oxygen consumption. The model is based on physiologic principles and data collected on human subjects. Properties of the model are discussed and its behavior is compared to that of human subjects under a variety of test inputs. A method of system identification is then developed. This method uses least squares and its generalizations to estimate parameters in NARMAX and Volterra series models. A useful set of input test signals is determined and a method for matching the properties of these signals to the system under test is discussed. The identification method is first applied to a Wiener system consisting of a first order linear subsystem cascaded with a "clipper" nonlinearity and the resulting NARMAX and Volterra series models are discussed. The identification method is then applied to the simulation model and NARMAX and Volterra models are constructed. These models are compared and their agreement with the original simulation is evaluated. The form and content of this abstract are approved. I recommend its publication. Signed Wi am Briggs
PAGE 4
CONTENTS INTRODUCTION 1 1 THE MODEL 6 The System to Be Modeled 6 The Computer 6 The Bicycle .. 8 The Measurement Devices 13 The Subject ......... 16 The Muscle Subsystem 18 The Venous Subsystem 34 The Lung Subsystem 38 The Cardiac Subsystem 39 The Visceral Subsystem 41 The Arterial Tree 41 The Overall System 42 An Example 43 Conclusion ...... ........ 50
PAGE 5
2 ANALYSIS OF THE MODEL 53 Introduction 53 Ramp Inputs 53 Step Inputs 59 Sine Wave Input 62 Constant Switching Pace Inputs 73 The Model with Noise 73 Summary ....... 80 3 DEVELOPMENT OF THE IDENTIFICATION METHOD 81. Introduction . 81 The Data Record 83 Candidate Models 93 The Volterra Se_ries 94 The NARMAX Model 98 The Assessment Rule 100 An Example . . 103 NARMAX Applied to the Example. 130 Volterra Series Applied to the Example 136 Conclusion . . . . . . . . . 139 4 IDENTIFICATION OF THE MODEL 142 Introduction . . . . . 142 v
PAGE 6
The NARMAX Model 146 The Volterra Model 172 Building Models in the Presence of Noise 180 5 SUMMARY AND FUTURE WORK 191 Summary of the Thesis. 191 Future Work 0 193 APPENDIX A 196 Wiener vs Hammerstein Models 196 BIBLIOGRAPHY 222 vi
PAGE 7
ACKNOWLEDGEMENTS The author would like to thank Carol Welch of the Denver Veterans Administration Hospital Pulmonary Function Laboratory and Charles Irvin of the National Jewish Hospital Pulmonary Function Laboratory who allowed us to use exercise testing equipment. Connie Borst and Paul Hogan, who work in the Veterans Administration Hospital Lab, and Steve Senn, Martin Carlos and Dave Gurka, who work in the National Jewish Hospital Lab have been very helpful. The author would also like to thank George Swanson, formerly of the University of Colorado Health Sciences Center who also allowed us use of exercise testing equipment. The author thanks his thesis committee and acknowledges the University of Colorado at Denver College of Liberal Arts and Sciences, which supported some of this work through a research assistantship. Vll
PAGE 8
INTRODUCTION The determination of V02maxthe maximum rate at which a subject can consume oxygenby graded exercise protocols has been in use for subjects rang ing from elite athletes to victims of severe cardiopulmonary disability for a long time [24]. In these graded protocols, the subject is presented with a gradually increasing work load and the rate of oxygen consumption is monitored until ei ther fatigue sets in, the subject can't maintain the load or the rate of oxygen consumption plateaus. This maximum rate of oxygen consumption (V02max) is then used both as a measure of conditioning state and to prescribe exercise regtmens. In recent years two possibilities have led to further interest in this system. One is the possibility of determining V02max by submaximal tests. This is par ticularly important in the post myocardial infarction patient or the patient with severe cardiopulmonary limitations where one would prefer to avoid maximal exercise. It is also important in limited, unusual environmentssuch as a space craft before reentrywhere assessment of conditioning is important but, exercise to maximal levels is inconvenient. The second possibility is that further insight
PAGE 9
into both the conditioning state and the underlying physiologic system may be obtained by studying the dynamics of the system under inputs other than gradual increase of the work load. In this thesis we will consider a broader problem which includes all of the abovethat of system identification for this input( work load )output( oxygen consumption) system. The specific system that we will consider is a human riding an exercise bicycle. A similar analysis would apply to a treadmill system or an arm cycle or any number of braked ergometer systems that exercise large muscle groups. We will use the bicycle system because there is a large installed base of test equipment available in hospitals and our analysis will be applicable not only to esoteric, expensive equipment in physiology labs, but also to equipment in hospital cardiopulmonary evaluation units. (There is also a large number of treadmills available, but their dynamic response is generally slower than the bicycles and, thus, they are not as useful.) The specific problem we consider is that of determining, from a given input time series of work load verses time and the corresponding output of vo2 verses time, how one can best construct a model that represents the system. We take the view that encoded in this mapping from the space of possible inputs to the space of possible outputs is all the information about the system that we can know. Once this system identification has been done, determination of V02max and the study of system dynamics is a relatively easy subproblem. More specifically we will limit ourselve to building models where the output of the system is expressed 2
PAGE 10
as a polynomial whose terms are products of previous system inputs and outputs; socalled lagged polynomials. As an example let the input be the time series x( t) and the corresponding output be y(t). A polynomial model of the type we are considering would be y(t) = ax(t1) + by(t2) + cy(t l)x(t 2). Such models have a long history in system identification [22]. If we restrict the polynomials to have factors involving the input only then the method becomes the impulse response method [25] in the linear system case and the Volterra method [39] in the nonlinear system case. Similarly, if we allow the lagged polynomials to include both inputs and previous outputs the method is the ARM AX model for AutoRegressive Moving Average with eXogenous inputs [11] in the linear system case and the NARMAX model for Nonlinear Autoregressive Moving Average with eXogenous inputs [27] in the nonlinear case. We will develop methods for deter mining the lagged polynomial model from inputoutput data. This will involve determining appropriate inputs with which to test the system and determining algorithms to construct the model from the resulting inputoutput pair. The approach we will use is a combination of analysis and numerical ex periment applied not to the original system, but to a simulation model of the system that we will construct working in MAT LABTA/.1 By using this technique it might appear that the result will be only as good as the model. We offer three 1 MAT LAB is a trademark of The Math Works, Inc. 3
PAGE 11
responses to this: The first is that the simulation model matches the response of the original system quite well, as we will show. The second is that the experimental design questions we will addresssampling time, input shape, etc.are the same whether using a model or the original system and so the techniques developed will be applicable even if the parameters or structures in the original system are not exactly those of the simulation. The third is that even if the simulation had no relation to any physiologic system it is still a complex, welldefined system and its identification is a fascinating problem in its own right that provides insight into both system behavior and identification methods. On the other hand, there are reasons for using a simulation instead of the original system. By using the simulation model we have been able to perform a large quantity of experimentation that would have been impossible with human subjects. We have also been able to conduct experiments such as impulse respouse and sustained high work loads that could not have been done with human subjects. By using the simulation model we have access to observablessuch as "true" transfer functions, venous pool sizes and instantaneous oxygen concentrations, and othersthat would not be available in a human experiment. Finally, the experiments are repeatable. While human performance varies from experiment to experiment and day to day, the simulation produces consistent results. The thesis is divided into five chapters. In the first we will develop the simulation model. In the second we will evaluate the simulation model and present 4
PAGE 12
qualitative and quantitative observations about the simulation and relate them to the original human system. In the third chapter we will develop methods of system identification and apply them to a simple system. In Chapter Four the methods of Chapter Three are used to identify the simulated system. Finally, in Chapter Five we will sumarize the results and present areas for future research. 5
PAGE 13
CHAPTER 1 THE MODEL The System to Be Modeled In this chapter we will build the simulation model. The physiological system that we are modeling is a human riding an exercise bicycle. The input to this physiological system is a work load determined by the braking force at the pedals of the exercise bike and the output is the oxygen consumption of the subject determined by measurement at the mouth. The physiological system, however, is only a part of the entire system we must model. The major system blocks are the bicycle, the subject and the measuring devices as shown in Figure 1.1. The experimental protocol is entered into a computer which then executes the experiment and controls data collection. We will include the computer in our analysis and so the full system is as shown in Figure 1.2. The Computer The data flow begins when the investigator either enters a protocol into the computer or chooses an established protocol already stored. This protocol consists of levels of work load and times to switch from one level to the next. In some available systems the number of switching times and levels is limited only by 6
PAGE 14
put Measuring Bicycle Subject Devices Out In put Figure 1.1: Basic block diagram of experimental setup. Input Protocol c D/A Bicycle 0 m p I u t I Measuring I Devices e r Output Figure 1.2: Block diagram of full experimental setup. 7
PAGE 15
computer memory and the time between switches is limited only by the computer cycle time. However, in sonie commercial units these choices are limited by the software sold with the package. We will explore the effect of the number of switches and the time between switches on the quality of identification in Chapter Three. The computer also administers the experiment in the sense that when the operator presses the start button and the experiment begins, the computer clocks the digitally stored protocol to a D I A (digital to analog) converter which sends an analog braking signal to the bicycle. It also begins to sample the ports of an AID (analog to digital) converter that interfaces to the output data lines for gas flow, gas concentration, heart rate and other variables of interest. The computer does the necessary numerical processing such as data smoothing and the calculation of oxygen volume from the measured values of total gas flow and percentage of oxygen in the total gas. It also provides information so the operator can follow the progress of the experiment. And finally, the computer stores the data values for further analysis and provides graphical output, either online or in later offline analyses. The Bicycle The bicycle receives the analog braking signal from the D I A converter. It has its own mechanical and electrical properties (for example, inertia) and acts as a subsystem with a known transfer function. The output of the bike subsystem is presented to the subject subsystem, as shown in Figure 1.2, as a force against 8
PAGE 16
which the subject must work. The subject maintains a near constant pedaling rate and so this force translates into a work load with units of watts. (We will use the MKS system of units.) We will model the bicycle as a linear system with known time constant T& (the time required for the system to reach 1 1/ e of its final value during a step input.) This must be determined by measurement and varies widely from bike to bike. In a high quality lab bike, T& may be as low as 1 second while in some commercial systems it is as high as 10 seconds. The equation governing the bike is assumed to be: (1.1) where y(t) is the load presented to the subject and i(t) is the voltage presented to the bike to determine braking force. Integrating ( 1.1) gives: y(t) = y(t0)exp((tt0)/r&) + (1/r&) {' i(s)exp((t s)/rb)ds. (1.2) ),0 In the simulation model there is an internal update time (or time step) that is chosen to be small relative to all time constants of the system being modeled. The switching times in the protocol are required to be multiples of this internal time. Thus, during any of these time steps the input to the bike is constantij during the j'h time step. The bike output for tin the ja' time step then becomes where tj is the starting time of the j"' time step. 9
PAGE 17
70 65 60 55 50 45 40 35 30\ I I l i \ .. r' '' :: .. . . .. / TIME(Secs) \ Figure 1.3: Input (solid line) and output (dashed line) of the bicycle. Figure 1.3 shows the output of the bicycle (dashed line) corresponding to an input (solid line) switching between 25 watts and 75 watts in a manner called pseudorandom, a commonly used test input for identification of linear systems. The internal time step for this run is 0.1 seconds and the time constant of the bike is 5 seconds. The dashed line in the figure represents the bike output and follows the "rapidthenslow" change expected from a first order linear system. There is another system component that should be discussed with the bicyclethe antialiasing filter. Since the data processing methods developed later will be digital, we are dealing with a sampled data system and thus the data must be filtered to avoid aliasing. (See [21] for a discussion of this topic.) One way to see the importance of this antialiasing filter in the setting of system 10
PAGE 18
identification is as follows: Suppose we input a sine wave whose frequency is below the Nyquist frequency1 into an arbitrary system and observe its effect. We have, in a sense, identified the system at that frequency. Now add the first harmonic of the sine wave to the input and suppose that it lies above the Nyquist frequency. Without the system will produce its response to the combined signal. However, when the data are sampled, the energy at the higher frequency will be folded onto a lower frequency and the identification algorithm will mistakenly interpret the resulting output as having come from the single sine wave below the Nyquist frequency. Clearly this is an unsatisfactory result and so we must filter the input as a function of the sampling frequency that will be used in the final identification algorithm. On the other hand, it is known [29] that an input signal that is used in an identification experiment must have "enough" frequency com ponents and so we must be careful that we do not remove too many frequencies with our filter. We implement this filter in the simulation as a fourth order Butterworth filter [42] placed early in the data path. Its effect is illustrated in Figure 1.4 which shows the effect of a varying sampling interval on an input that switches every ten seconds between 25 watts and 75 watts. As one goes down the figure the sampling time goes from 0.5 second to 1 second to 2 seconds and finally to 5 seconds. One can literally see the high frequencies disappear from the signal. Although the effect on the output of the system depends on the frequency 1The Nyquist frequency is the lowest frequency such that there are two samples per cycle and, thus, its value in Hertz is one half the sampling rate in samples per second. 11
PAGE 19
00 !E 00 !:: " 100 50 0 0 100 50 0 0 100 50 50 50 50 tOO !50 200 250 TIME(SECS) !00 150 200 250 TIME(SECS) J !00 !50 200 250 TIME(SECS) 100,. 50 0 50 100 !50 200 250 TlME(SECS) Figure 1.4: Effect of sampling and filtering on an input. The input is a step shown in the top plate. In the successive plates the sampling time increases and the corresponding antialiasing filter becomes more narrow. 12
PAGE 20
i 1 response of the components, it is clear, as we will discuss later, that this effect is important. The actual realization of the filter in the true system can be in hardware between the computer and the bicycle or it can be realized digitally in the computer. Causality is not required in these experiments because inputs are determined entirely before the experiment begins. The way we will implement this filter in the simulation is with the wellknown signal processing method of passing the input through the Butterworth filter and then reversing the output and passing it again through the filter. (This is done for ease of implementation. If we have a hardware filter or any filter with nonzero delay we could easily correct for it.) The Measurement Devices Returning to the block diagram of Figure 1.2 we will discuss the Measurement Devices block. In a typical experiment the raw variables that are measured are instantaneous heart rate, respiratory rate, gas flow at the mouth, oxygen consumption and carbon dioxide concentration at the mouth. Many systems provide for simultaneous measurement of other physiologic observables such as hemoglobin saturation by ear or finger oximetry. The heart rate can be measured by counting the time between voltage peaks of the QRS complex of the EKG and inverting this to give beat to beat heart rate. In the systems we have used this value is available as a voltage output from the EKG monitor and is input to the computer on one of the channels of the 13
PAGE 21
A/D converter. The computer samples the channel and may or may not smooth it before storing it depending on the experimental protocol. Gas flow is typically measured by a fan of low moment of inertia placed in the gas flow path just distal to the mouth. The fan rotations are counted by a photocell sensor that responds to reflected light pulses. The output current of the photocell is used to charge a capacitor that acts as an integrator. The capacitor is discharged with each reversal of the fan and so the instantaneous voltage across it corresponds to the volume of gas passing the fan so far in each breath. This voltage across the capacitor is also input to a channel on the A/D converter and the computer samples it according to the experimental protocol. By measuring the time between reversals, the computer can determine instantaneous respiratory rate and by differencing this total volume signal the computer can determine the volume of gas flow per sampling interval. Some lab systems use a mass spectrometer to determine the composition of the expired gas. These systems output analog signals corresponding to percentage oxygen, percentage carbon dioxide and percentage nitrogen to channels on the A/D converter. The percentage oxygen signal is sampled and multiplied by the corresponding gas flow signal in the computer to yield oxygen flow. This signal is integrated over an exhalation and compared to the value for the corresponding inhalation to get the breathtobreath oxygen consumption. A similar calculation may be done with the carbon dioxide data to get the breathtobreath carbon dioxide production. 14
PAGE 22
I Other systems perform carbon dioxide analysis with an infrared absorption analyzer and oxygen analysis with a fuel cell. It is important to calibrate these measurement systems and to determine their inherent time delays. However, once this is done we believe that they accurately reflect the physiologic variables and we will model them in the simulation as a straight passthrough from the state variables output by the Subject Box in Figure 1.2. VIle must mention, however, that the output from the system is selfsampled on a breathtobreath basis. Although we have access to flow data, gas concentra tion and other variables at arbitrarily high sampling rates (some systems sample these in the kilohertz range), these data are used only to calculate how much oxy gen is inhaled and how much is exhaled as discussed above. These are subtracted to get one number per breath that represents the oxygen consumption during that breath. We have found that a subject can pace respirations to a nearly fixed rate by varying only the depth of respiration within limits. There has also been work done regarding best breathtobreath interpolation techniques. There are many implementations of these methods and a real bike experimental system can provide output time series sampled on a breath to breath basis or, by interpo lation, at any sampling rate. Our simulation model will produce observations of these variables as frequently as one observation per time step. 15
PAGE 23
The Subject Now consider the Subject Box in Figure 1.2. The input to this box is the work load presented by the braking force applied to the pedals and the output is the vector of sampled physiologic observables as discussed in the previous sec tion. Figure 1.5 shows the major anatomical components of the cardiorespiratory system: L represents the lung and upper airways that are the source of external oxygen for the system and the sink for carbon dioxide; C is the heart that pumps the blood that carries the oxygen to the tissues; M represents the muscle compartment; 0 represents the "other" consumers of oxygen (brain, viscera, etc.); V represents the venous system that returns the oxygen depleted blood from the muscles and viscera to the heart and lungs; and A is the arterial system that carries the oxygenated blood from the heart and lungs to the peripheral com partments. In this diagram one can trace the oxygen path from external air to mitochondria that will eventually determine V02 (the rate of oxygen consump tion). U(t) is the work load input which will act on the muscles and the cardiac subsystems as represented by the arrows in the figure. Y(t) is the observed out put vector that acts as input to the Measuring Devices Box. We will now discuss and model each subsystem in turn. 16
PAGE 24
/Y(t) L A v M l:J( t) Figure 1.5: Major anatomical components of the cardiorespiratory system. 17
PAGE 25
MUSCLE BLOOD FLOW WORK LOAO VENOUS OXY GEN INPUT CONCENTRATI ON MUSCLE OXYGEN L CONSUMPTION 1MUSCLE SUBSYSTEM Figure 1.6: Overall model of the muscle compartment. The Muscle Subsystem For our purposes U(t) (exercise load) will affect the muscles in two ways. It will cause oxygen to be locally consumed as glycolysis occurs to fuel the muscle contractions. It will also regulate local blood flow to carry oxygen and other nutrients into the muscle and to carry carbon dioxide, lactate and other metabolites away. We will model these as two independent processesone with work load as input and muscle blood flow as output and a second with work load as input and muscle oxygen consumption as output as in Figure 1.6. This type of model has been used before. ([5], [14].) A fundamental property of linear systems is superposition which implies that doubling input doubles output. It is clear that the two processes discussed above are nonlinear for large input variations, since if we repeatedly double the work load there will be a load 18
PAGE 26
beyond which the flow no longer doubles. While the image is not so dramatic, the local consumption of oxygen will also stop doubling as the load doubles beyond some limit. It is also clear that these subsystems have dynamics in that the response to an input is not instantaneousone does not note shortness of breath or increased heart rate immediately after increasing work load. Much research has been directed toward determining the dynamics of these subsystems as dis cussed in [5]. We will assume that they are first order. The simulation model is modular and the implementation of higher order kinetics would be straight for ward. One method for modeling nonlinear systems with dynamics is to isolate the nonlinearity in one subsystem with no dynamics(ie. instantaneous response and no memory) either preceded or followed by a second subsystem that models the dynamics and is linear [22]. If the linear subsystem precedes the nonlinear static subsystem the model is called a Wiener model and if the nonlinear subsystem precedes the linear subsystem the model is called a Hammerstein model as shown in Figure 1. 7. We will represent both of the muscle level subsystems (work load to blood flow and work load to oxygen consumption) in one or the other of these ways. In the case of the blood flow subsystem, for the Wiener model the work load will be input to a first order linear system. Since the simulation model is discrete in time we will conceptually sample the input work load at each time step and the model will evolve between consecutive internal steps governed by the first order linear model. The output of this subsystem will then become input for the nonlinearity 19 _LL ________ __
PAGE 27
LINEAR SUBSYSTEM I NONLINEAR SUBSYSTEM WIENER MODEL NONLINEAR SUBSYSTEM 1,_1 1 L_ ____________ LINEAR SUBSYSTEM HAMMERSTEIN MODEL Figure 1.7: Wiener and Hammerstein models. to yield the overall output. On a physiologic level, this means that the input work load is delayed and its power density spectrum is modified by some underlying mechanism. However, this is quite reasonable since we are seeing the modulation of local vascular resistance presumably mediated by local messengerseither chemical or neurologicalthat do not respond immediately or uniformly to work load. This "apparent force" that is trying to alter flow is then modulated by the reality that the arteriolar j capillary system has finite capacityas represented by the horizontal asymptote on the function representing the nonlinear portion of the model. This is seen in Figure 1.8 which shows such a Wiener system. In the Hammerstein model the cascade is reversed and the argument is that the input work load is modulated by the nonlinearity first to yield an "apparent force" which becomes the input to the linear dynamic subsystem. Similar heuristics apply to the oxygen consumption subsystem. In this case the forces are not the opening of blood vessels that lead to changes in flow, but 20
PAGE 28
J LS I I L_ L_ _____ _J Figure 1.8: Wiener model with limiting nonlinearity. LS is the linear system and the saturating curve represents the nonlinear system. are derived from chemical potentials that drive reactions down redox chains and lead to changes in oxygen consumption. We must now decide which model to useWiener or Hammerstein. Figure 1.9 shows the result of a sine wave input with high frequency relative to the time constants of the system to the Wiener model and Figure 1.10 shows the result when the same sinusoid is input into the Hammerstein model. The nonlinearity in each is the saturating type shown in Figure 1.8. The output in each case is nearly a sine wave with the same frequency as the input. However, the mean of the output for the Hammerstein model is less than for the Wiener model. Figure 1.11 shows a similar plot for a sine wave input with low frequency relative to the time constant for the Wiener system and Figure 1.12 is the same for the Hammerstein system. Note that here the outputs from both models are periodic and of similar appearance, but nonsinusoidal. Again the mean for the Hammerstein case is less than for the Wiener case. Thus, although the models produce outputs that are 21
PAGE 29
43.5 43 1\ 42.5 z ::;, f... 42 ::;, 41.5 \i 41 40.5 H II II I I I I I I I I I I I I I I I ' I I I I /ill \1 I I I Ji I II I I \ I 0 50 100 !50 200 250 300 TIME(SECS) Figure 1.9: Output for high frequency sinusoidWiener model. 22
PAGE 30
44 42' I Ul r .... z ::0 I .... 41 I H ::0 0.. .... ::0 I 0 10 I I 39 1 38 0 50 100 150 200 250 300 TIME(SECS) Figure 1.10: Output for high frequency sinusoidHammerstein model. 23
PAGE 31
46 42 nn 1 "' I I \ I \ J ... z :::> \ I I \ ... 40 \ I :::> I 0.. \ I ... :::> 0 38 v :J 0 50 100 150 200 250 300 TIME(SECS) Figure 1.11: Output for low frequency sinusoidWiener model. 24
PAGE 32
46 44 42 \ I I UJ \ \ l ,.. z 40 "' ,.. "' 0.. 38 ,.. 0 I 0 I 36 34 32 0 50 100 150 200 250 300 TIME(SECS) Figure 1.12: Output for low frequency sinusoidHammerstein modeL 25
PAGE 33
similar, the levels are different and the nonlinear and linear subsystems do not commute implying that the two models are not equivalent. An analysis of this phenomenon is instructive, but it does not resolve the choice between the two models and so is presented in Appendix A. The model we will use for both subsystems of the muscle system is the Wiener model. In fact, we can get either model to match human subject data by adjusting the inputoutput curves associated with the nonlinearity. However, we feel that the linear system forced by the work load followed by the nonlinear "capacity limiting" function represents the physiologic reality better than the reverse order of the Hammerstein model and so we will choose the Wiener model. The linear subsystems will all be simulated as discussed in The Bicycle section. However, here we do not assume that the input remains constant during the internal time step and we approximate the integral in the evolution equation ( 1.3) with the trapezoidal rule to get: y(ti+l) = y(tj)exp( (flt/r;)) + (1/r;)(f1t/2)(U(tj)exp( (LH/r;)) + U(ti+I)) (1.4) where y is the output of the linear subsystem, f1t = ti+l tj, [tj, tj+J] is the j'h internal time step and T; is the time constant for the respective subsystem. (In the following Tqm will be the time constant for the muscle flow subsystem and Tvm will be the time constant for the muscle oxygen consumption subsystem.) There are some clinical data regarding the actual values of Tqm and Tvm ob26
PAGE 34
tained from isolated preparations of animal muscles and indirect observations in man. There is also no a priori reason to assume that the "off" kinetics (response to decreasing work load) are the same as the "on" kinetics. (See [23].) Barstow and Mole [5] construct a model to represent V02 under the influence of step inputs. They use values of Tqm from one to forty seconds and values of Tvm from ten to forty seconds and we will also consider values in this range. We will also allow the "on" and the "off" time constants to be different. It is also possible that there are actual delays in the system. (See [23].) These are modeled in the simulation by simply shifting the time at which the input is applied. Of course, all of these constants vary from one experimental subject to anotherjust as each experimental subject has their unique set of parameters, so the simulation corresponding to that subject will also have its unique set of parameters. The next question to be answered is how to model the zeromemory, nonlinear compartments. Note that when the input is constant the linear compartment output approaches a steady state value as does the output of the cascaded model. There are data available relating this steady state input to the steady state output [31] and we will use these results to construct curves of both oxygen consumption and muscle flow verses input load for a general subject. We will then use these curves as the nonlinear part of each subsystem. Again, these curves are unique to each individual and must be determined for each experimental subject we will do in Chapter Two. However, Figure 1.13 is a plot of the relation between muscle oxygen consumption and apparent load presented by the output of the 27
PAGE 35
50 45 0 40 0 "" 0 35 :g 30 0 :;:l "" 25 s 20 0 0 N 15 0 0 10 0 ..., 5 0 0 50 100 150 200 250 300 350 400 450 500 Work load (watts) Figure 1.13: Nonlinear relation between work load and muscle oxygen consumption. linear dynamic subsystem constructed for a hypothetical subject from data in several sources such as [31] and [10]. Figure 1.14 is a similar plot of muscle blood flow verses apparent load. In each there are portions that are nearly linear, especially if we limit our considerations to small changes in the work load, and each approaches a horizontal asymptote with increasing work load, as expected. Note that we do not know the exact shape of these curvesthey are constructed by interpolating and extrapolating data. Their shape suggests that they may fit 28
PAGE 36
I I "' '! 0 r;:: 250 ; 200 I 150 100 50 100 150 200 250 300 350 400 450 500 Work load {watts) Figure 1.14: Nonlinear relation between work load and muscle blood flow. 29
PAGE 37
a MichaelisMenton curve: k1 *input output= k3 + (k )" 2 + tnput There is some physiologic support for this in that we can view either subsystem as flow through a limiting compartmentspatial limiting in the flow case and limits on the rate at which oxygen can be consumed in the oxygen case. Similarly, the MichaelisMenton relation attempts to model reaction kinetics as flow through a system with limited enzyme available. At any rate, we are able to fit the empirical data with the MichaelisMenton model and we will use them in the simulation model. It can be argued that the true purpose of the identification algorithm is to determine the shape of these flow verses load and oxygen consumption verses load curves for all load values. Suppose, for example, that we want to determine V02max from submaximal inputs. Suppose further that the rate limiting factor is flow. Then the input must be designed to determine, at least implicitly, how flow extrapolates to its maximum. If the true shape of flow verses load in the nonlinear subsystem were to be as in Figure 1.15, even the best designed test at loads below 80 watts would give little hope of determining V02max Thus, in summary, we have modeled each muscle level subsystem as a Wiener model. Figure 1.16 shows the result of a piecewise constant input into such a. model. The top plot shows the input of work loa.d verses time. The second plot shows how the linear system responds to this a.nd finally the third plot shows 30
PAGE 38
0 u T p u T ANOTHER POSSIBLE INPUTOUTPUT RELATION 80 INPUT(WATTS) Figure 1.15: Another possible nonlinear relation between work load and muscle blood flow. 31
PAGE 39
INPUT TO LINEAR SYSTEM I II_ L I + Yi si FINAL SYSTEM OUTPUT 1( nu I _, OUTPUT OF LINEAR SYSTEM NONLINEAR SYSTEM Figure 1.16: Overall summary of Wiener model. 32
PAGE 40
how the nonlinear system attenuates the output of the linear system to produce the overall system output. This last plate shows both the output of the linear system and the final system output superimposed. There is one last step to complete the modeling of the muscle system. Recall from Figure 1.5 that the Muscle System also has input from the arterial tree and produces output into the venous system. The state variables that we will observe in both the arterial and venous sides is oxygen concentration in milliliters of oxygen per milliliter of blood. Since we really want milliliters of oxygen to be a measure of oxygen quantity, we could use a mass measurement or even count molecules, but the tradition is to correct the volume of oxygen per unit of blood to standard temperature and pressure and call the units milliliters percent. From the above discussion we regard the local muscle blood flow and the local muscle oxygen consumption as a function of exercise load level. From this we can determine the relation between oxygen concentration on the arterial inflow side and oxygen concentration on the venous outflow side by a statement of conservation of mass called Fick's principle. (See [40].) In this case it is: (Caoz(t)Cvoz(t)) Q(t) = Qcn(t) (1.5) where Q is the local flow (ml/sec), Caoz is the arterial concentration of oxygen (ml oxygen/ml blood), Cvo2 is the venous oxygen concentration (ml oxygen/ml blood) and Qoz is the local oxygen consumption (ml oxygen/sec). It is known that an exercising subject tries to maintain several arterial values m a narrow 33
PAGE 41
range as long as the inspired gas mixture is not varied [40]. One such value that is defended is arterial oxygen content at about 21 ml/ 100 ml blood [5]. We have found by ear oximetry that this value does in fact decrease during vigorous exercise. However, we will assume it to be a constant in the simulation and so we can solve equation (1.5) for Cvo2 as a function of local flow and local oxygen consumption and, thus, of input work load. This then becomes an input to the Venous System, which we model next. The Venous Subsystem In view of the previous discussion, it is clear that as exercise load changes two things will occur: Local flow into the venous system will change and the blood that flows into the venous system will have varying oxygen content. We will view the venous system as in Figure 1.17. There are two sources of inputone is flow and oxygen content from the muscle box and the other is flow and oxygen content from the visceral (or "other" compartment). (Recall Figure 1.5.) The output of the venous system is the flow and blood oxygen concentration into the right side of the heart. It is well known that the venous system can act as a variable size storage reservoir and we will model this as three fixed volume pools and one variable volume mixing chamber. Two of the fixed size pools connect the source at the muscle box and the source at the visceral box with the mixing chamber and the third fixed size pool connects the mixing chamber with the sink at the right heart. We will we not take into account possible beat to beat blood storage 34
PAGE 42
OUTPUT TO RIGHT HEART=====l """ "" eoo" 17 INPUT FROM ::::::====1 OTHER SYSTEMS VARIABLE SIZE MIXING CHAMBER FIXED SIZE POOL ,INPUT FROM _____ __j MUSCLE SYSTEM __.::.:__ _______ Figure 1.17: The model for the venous subsystem. 35 ELEMENT OF BLOOD
PAGE 43
m the pulmonary vasculature and will require this flow into the right heart to equal the left heart cardiac output. The size of the venous pool has been determined by several methods and we will use standard values for the equilibrium values and apportion them among the four in a manner that gives good agreement with experiment. The input (flow and oxygen concentration) from the muscle compartment is a known function of work load as derived in the previous section. The flow and oxygen concentration from the visceral compartment will be determined by the methods presented in the subsection titled The Visceral Compartment. Some previous models of the cardiorespiratory system have assumed that blood from the muscle compartment is fully and instantaneously mixed with the venous pool. However, we will track small elements as they exit the muscle compartment and traverse the venous pool to the mixing chamber, treating this part of the venous system as a first infirst out buffer. This part of the system between muscle and the mixing chamber will be considered a constant volume pipe and each discrete bolus of blood will have associated with it a content of oxygen that remains constant as the element traverses this part of the venous system. The pipe will empty into the mixing chamber at a rate determined by the flow into the pipe, which is simply the output of the muscle compartment. The input from the visceral compartment will be treated in a similar manner. I In the mixing chamber we will assume complete mixing occurs. Since we know the initial size of the chamber and since we know the flow into it and flow out 36 j
PAGE 44
of it (cardiac output) we can determine the size of the mixing chamber at any time and it will expand and contract depending on the various time constants of the flow systemscardiac output, muscle flow and visceral flow. Also we will assume that the initial oxygen concentration in the chamber is the average rest oxygen concentration of venous blood, .16 ml oxygen per ml blood. (See [31].) Thus, since we know the concentration of oxygen flowing into the chamber we can calculate the total oxygen content of the chamber and the oxygen concentration. This oxygen concentration is then associated with each element of blood in the mixing chamber at the current time. These elements then become the output of the mixing chamber and we view them as being pumped out of the chamber each time step in a quantity equal to the cardiac output. We then track these elements as they traverse the fixed volume pipe leading from the mixing chamber to the right side of the heart and then into the lung and then into the left side of the heart for pumping into the arterial tree. Note that the transit time is flowdependent in that the time delay (TD) in each part of the system is determined by an integral of the flow: {' flow dt = volume ltTD In the simulation we model this by tracking elements of blood through the venous system as time evolves discretely. A continuous form of this model for the venous system appears in the literature as early as the 1967 reference by Grodens [20] and in a more recent model by Essfeld and Hoffman [14]. We feel it is much more 37
PAGE 45
realistic than the complete mixing modeL It also introduces a variable delay nonlinearity into the system that is fascinating in its own right. The Lung Subsystem As we discussed briefly in the previous section we do not allow variable storage of blood in the pulmonary vasculature. During each time step of the simulation a quantity of blood enters and the same quantity of blood exits the lungs and both of these equal the cardiac output. The input to this subsystem comes from the fixed pipe connecting the venous mixing chamber with the right heart as shown in Figure 1.17 and the output of this subsystem is the oxygenated blood that is input into the left heart. This subsystem calculates the simulation output variable that we are most interested inV02 Again we use Fick's principle and make the assumption that the arterial oxygen concentration has a fixed value of 21 ml oxygen per 100 ml of blood. In its appropriate form Fick's principle is, using the previous notation: or, since arterial oxygen concentration = .21 ml oxygen/ml blood, V02 = (.21 Cvo2) Q Since we know the venous oxygen concentration and cardiac output we can calculate 1f02 and we finally have this value as a function of work load input. 38
PAGE 46
It is now possible to determine ventilation or other easily observable respiratory parameters by making assumptions about lung mechanics. We will not do this since we are interested in V02 However, we must note that although we have oxygen consumption measurements available at time samples limited only by the time discretization in the simulation, in reality one is practically limited to observing the process only on a breathtobreath basis and the observation is done at the mouth not at the alveoli. We will assume full mixing in the lungs and assume a constant dead space. This makes our alveolar values the same as the values collected at the mouth. We resolve the breathtobreath question either by integrating our nearly instantaneous observations over a breath cycle to produce breath to breath output or by assuming there is an interpolation algo rithm, as discussed earlier, that essentially turns breathtobreath sampled data into continuous data. The Cardiac Subsystem Cardiac output will be modeled in a manner similar to muscle blood flow. We again use the Wiener model of a linear subsystem with first order dynamics cascaded into a zero memory nonlinear subsystem. The input to the system is work load and the output is the cardiac output, which as discussed before, is the amount of blood pumped from the venous system through the lung system and into the arterial system on each internal time interval of the simulation. More complex models for the linear subsystem have been suggested (see Fugihara [16]) 39
PAGE 47
400 350 300 Q 250 :_ 0 200 u 150 100 50 0 50 100 150 200 250 300 350 400 450 500 Work load (watts) Figure 1.18: Nonlinear function between work load and cardiac output. and since the simulation is modular they would be easy to implement. In this thesis we will use only the first order system with time constant Teo in the one to forty second range and no fixed delay. Like the time constants for the subsystems discussed earlier the values of Teo vary from subject to subject. The nonlinear subsystem is similar to the nonlinear subsystem for muscle flow. The known steady state values from [31] or [10] are used to generate a table of cardiac outputs as a function of work load. As before, the intermediate values are determined either by linear interpolation or by fitting to a MichaelisMenton equation. (See Figure 1.18.) 40
PAGE 48
The Visceral Subsystem This part of the system is represented as 0 (for Other) in Figure 1.5. It represents such body parts as viscera, brain, nonexercising muscle, skin, and other physiologic compartments not otherwise represented in the model. We model it by saying that the flow into it is the cardiac output minus the flow to exercising muscle. We assume that its consumption of oxygen is relatively fixed and independent of work load. At rest this compartment uses about 300 ml oxygen per minute or 5 ml oxygen per second [31]. The oxygen content at the input to this compartment is the arterial value of .21 ml oxygen per ml blood by the same argument used at the input to the muscle compartment. We can, thus, apply Fick's principle again and determine the oxygen concentration at the output. We assume no storage in this system and so the flow at the output is the same as the flow at the input. This output flow and output oxygen concentration become the input to the fixed venous compartment connecting the visceral compartment to the mixing chamber. The Arterial Tree The final part of the system is the arterial tree that links the cardiac subsystem with the visceral compartment and the muscle compartment in Figure 1.5. The model we use is quite simple in that we will view it as a fixed volume pipe filled with blood at the constant oxygen concentration of 21 ml oxygen per 100 ml 41
PAGE 49
blood, as discussed earlier. This approximation allows us to essentially ignore time delays and flow on the arterial side. The Overall System Now that we have constructed each subblock of the overall system let us again consider the overall system as an inputoutput system. As we have seen, one possibility is to view the input as an element of R" where the j'h component of this vector is the value of the input during the j'h time interval. Each element is held for one switch time to become the input to the subsystem downstream. The inputoutput mapping in this case is from R" to the output space. Since we know the output of the computer (input to the bicycle) at all times it is possible to conceptually resample it at any rate that we choose including one that is not equal to the time separating the components of the original input vector. This is the view most commonly found in the physiology literature where, for example, a total time of experiment might be 10 minutes (600 seconds). The switching time could be taken to be 5 seconds (ie. the computer would switch levels every 5 seconds) and son= 600/5 = 120. (This is then in the previous paragraph.) The output of the sample and hold is then resampled at 1 Hz., thus yielding one sample every second and so m = 600/1 = 600. We can, in this case, view the experiment as a mapping from Rm to the output space. A third possibility is to view the original R" input as a sample from a finer grained, even continuous, "preinput". In this case the computer can be consid42
PAGE 50
ered a sampler in series with a zero order hold [15] and, in this case, the model mapping is from a subspace of a function space limited only by the restrictions imposed on its frequency content by the sampling theorem. The output can be viewed in a similar way. One can accept the breathtobreath values actually output from the physiologic system, in which case the output space is R\ where k is the number of samples taken of the output. Or one can smooth them as discussed in the section about the measuring devices. If this smoothed output is used directly, the output space is a function space. If the smoothed output is resampled, the output space is then R", where p is the number of resamples. These possibilities will become important in Chapter Three where we will apply methods to identify systems from inputoutput data and they will be discussed further there. Also, since it is really the subject subsystem in Figure 1.3 that we are interested in, and since we know the transfer functions for the bicycle subsystem and the measuring devices subsystem we can "invert" these subsystems and study just the mapping from the input to the subject subsystem to its output. Along these same lines we will consider the question of what input signal can be input to the overall system to produce a desired testing signal at the input to the subject subsystem. An Example In tills section we will present an example of how the simulation works. Since we feel that the model represents the physiology we will not add the word simu43
PAGE 51
!50 UJ .... z 100 ;:o 50 TIME(SECS) Figure 1.19: Work load in watts that is input to the example system. lated in the following and, for example, say "muscle blood flow" for both muscle blood flow in the human experimental system and the simulated muscle blood flow in the simulated system. Figure 1.19 shows a sample input choosen at random, subject only to the requirement that its amplitude range is limited to the physiologic range of zero watts to about 250 watts. In the real world experiment, the investigator decides how frequently to apply the sample and hold to the signal (say every ten seconds as we will use here) and he enters these corresponding values into the computer. When the experiment is started these values are clocked out to the bicycle through the D/ A converter. Figure 1.20 shows this output. As we discussed earlier we will process the data digitally and so we must pass the sampled and held data through an antialiasing filter. In this case we will sample at 1 Hz. and so the Nyquist frequency is 0.5 Hz. Figure 1.21 shows the result after the antialiasing filter. Since few high frequencies were present in the original signal and we are sampling at a relatively high frequency we see little effect from 44
PAGE 52
200 150 (fl ' t: J z 100 p ' _, ' 50 : j .J \ ' J 0 0 50 _! ' ' ' ' 100 ' , _: ' TIME(SECS) , : 150 ' ' \ : , i, \ _, ' :1 ', \ j 200 _, , ' 250 Figure 1.20: Original input (solid) and the sampled and held signal (dashed). 150 100 {\ 51 Vl( 0 50 100 150 200 250 TIME(SECS) Figure 1.21: Output of the sample and hold after it is passed through the antialiasing filter. the filter except a hint of the Gibbs phenomenon. We assume that the bicycle subsystem is linear with time constant rb = 1 sec. The output of the bike, which is the load the subject actually rides against is shown in Figure 1.22. This input acts on the three Wiener subsystems and produces the corresponding outputs. Figure 1.23 shows the resulting muscle blood flow verses time. If we compare this output with the original input we see a time lag as we would expect and we can also see the saturation effect of the nonlinearity. Figure 1.24 shows 45
PAGE 53
!50 100 50 0 50 100 !50 200 250 TIME(SECS) Figure 1.22: Output of the bicycle subsystem. 300,..,.. 250 lfJ t: 200 z ;:, !50 0 50 100 150 200 250 TIME(SECS) Figure 1.23: Muscle blood flow in ml. per second resulting from the input shown in Figure 1.19. 46
PAGE 54
40 Ul .... 30 z ::> 20 0 50 100 150 200 250 TIME(SECS) Figure 1.24: Muscle oxygen consumption in ml. oxygen per second resulting from the input of Figure 1.19. the resulting muscle oxygen consumption verses time. Note that the muscle blood flow and the muscle oxygen consumption curves are similar. This follows since the time constants for the two systems are chosen to be the same in this example and since the curves representing the nonlinearities are also similar. Figure 1.25 (top plate) shows the resulting cardiac output verses time plot from the input in Figure 1.19. Again it is similar to the two previous plots. The bottom part of Figure 1.25 shows the output of the linear subsystem of the Wiener subsystem for cardiac output. Comparing the overall output (top) with the intermediate output of the linear system we see the effect of the modulating nonlinearity. The muscle blood flow and muscle oxygen consumption outputs combine as discussed earlier to produce the oxygen concentration of blood entering the venous pool and a plot of this quantity is shown in the top part of Figure 1.26. The corresponding concentration of oxygen in the venous pool is shown in the bottom 47
PAGE 55
350 300 "' t:; 250 z ::> 200 150 0 50 100 150 200 250 TIME(SECS) 200 150 "' t:: 100 z ::> 50 0 0 50 100 150 200 250 TIME(SECS) Figure 1.25: Cardiac output resulting from the input in Figure 1.19 (top) and the corresponding output from the linear subsystem of the cardiac output Wiener system (bottom). 48
PAGE 56
0.1 (/) '"' J z ::0 0.08 0.06 0 50 100 150 200 250 TIME(SECS) 0.12 0.1 (/) '"' z ::0 0.08 0.06 0 50 100 150 200 250 TIME(SECS) Figure 1.26: Oxygen concentration of blood entering the venous pool (top) and the oxygen concentration of blood in the venous pool (bottom) resulting from the input in Figure 1.19. Both graphs have units of ml. oxygen per ml. blood. 49
PAGE 57
50 40 30 TIME(SECS) Figure 1.27: Oxygen consumption of the system in ml. oxygen per second resulting from the input work load in Figure 1.19. plot of Figure 1.26. Note that the result is smoother than the top plot and lags behind it, both a result of the fact that the mixing chamber is an integrating device. The blood is then pumped from the mixing chamber and past the lungs where it takes up oxygen. This determines the final oxygen consumption and, thus, the final output of the system and is plotted in Figure 1.27. Conclusion At this point we have developed a computer simulation model of the human cardiorespiratory system and an associated experimental environment that takes as input a collection of numbers that represent work load and an associated collection of times to switch from one work load to another. The model produces an output consisting of observations of 1f02 and several internal variablessuch as cardiac output, venous oxygen concentration, and venous volume. If we assume 50
PAGE 58
the dynamics are first order we have fifteen parameters. Twelve of these are represented as an "on" and an "off" time delay and an "on" and an "off" time constant for each of the three linear subsystems discussed above. We use "on" here to mean that the input to the system is increasing and "off" to mean that the input to the system is decreasing. The three other parameters are the sizes of the fixed venous compartments. We also have initial conditions for each of the systems and the initial size of the mixing chamber. The three curves of steady state values verses work load are also variables in the simulation. The system has nonlinearities that arise from several sources. One is the application of Fick's principle; for example, the calculation that determines venous oxygen content, requires the quotient of oxygen consumption and blood flow and so introduces a quotient nonlinearity. Another application of Fick 's principle is in the calculation of oxygen consumption at the lung. Here it requires the product of venous oxygen content and flow and so introduces a product nonlinearity. Other nonlinearities arise from the nonlinear portion of the Wiener model used to model cardiac output, muscle blood flow and muscle oxygen consumption. A third nonlinearity arises from the complicated, time delay that occurs in the venous system. We believe the idea of using a Wiener model with the nonlinearity represented as we have done here is new to cardiorespiratory modeling. It allows for complicated inputs and also allows different on and off dynamics. It also establishes that an important part of the identification of the system is the determination of 51
PAGE 59
the inputoutput curves for the nonlinear subsystems. We also believe that the venous model where each bolus of blood is tracked across the venous pool is also new. In the next chapter we will present and discuss sample experimental runs and apply some techniques from systems analysis to the simulation. In Chapter Four we will use the simulation to apply and compare time domain system identifica tion techniques and show how these results have application in the actual testing of the human subject. 52
PAGE 60
CHAPTER 2 ANALYSIS OF THE MODEL Introduction In this chapter we will determine the output of the simulation under several inputs often used in system analysisramps, steps, sine waves and a class of signals often called constant switching pace (CSP) [30] that are useful in testing and identifying systems. We will discuss why the simulation responds as it does and compare the simulated output with output from human subjects. In the last section of this chapter we will introduce a noise model into the simulation. Ramp Inputs The first input that we will consider is a ramp as shown in Figure 2.1. This is the type of input used clinically to determine V02max In such an experiment the braking voltage to the bike is slowly increased (and so the work load) until the V02 plateaus, the subject fatigues or stops for some other reason. Figure 2.2 shows an actual response from a human subject. In this plot we have averaged the breathbybreath output over eight consecutive breaths using a standard proto col used in exercise testing [24] and so this output appears relatively smooth. Figure 2.3 shows the human response with the simulation response superim53
PAGE 61
200 (/) "" ;:, 100 < 3 0 50 100 150 200 250 TIME(SECS) Figure 2.1: Sample ramp input of work load vs time. 2500 z S1 2000 ....._ ..., :::< {:;;' "' < 1500 .... "' ::> z "' "' 1000 ,.. >< 0 500 0 50 100 !50 200 250 LOAD(WATTS) Figure 2.2: Output of a human subject under the ramp input. 54
PAGE 62
50 u "' /' {/) ,'', "/..., /'' v' ;:;',/ ::::> / "" J ,.. / ::::> v 0 0 0 50 100 150 200 250 300 TIME(SECS) Figure 2.3: Human subject output (dashed) and simulation output (solid) under a ramp input. posed. They agree quite well and, in fact, we can make them agree very closely because we use this test to fit each of the cardiac outputwork load, muscle oxygen consumptionwork load and muscle blood flowwork load curves discussed in Chapter One to the given human subject. These three curves determine the vo2 at each input level and we could determine them by testing the subject with many different constant inputs. Instead, we choose the three curves above by trial and error in such a manner that they produce output that matches the ramp. Note that since there are three curves at each data point we have three degrees of freedomone for each curve. This process only uses up one degree of freedom at each point and so we have two left to fit further subtleties in the data if needed. One thing to note is that the output is nearly linear over a large input range. This is sometimes cited as evidence that the V02 system is linear. In reality, this 55
PAGE 63
apparent linearity is an artifact of the experiment. Recall from Fick's principle that the oxygen consumption at the lung is determined by the cardiac output and the oxygen concentration gradient between the blood flowing into the lung and the blood flowing out of the lung (recall Chapter One): VOz = (.21c\'02) *co, where .21 is the oxygen concentration of the arterial blood, C,02 is the oxygen concentration of the venous blood and CO is the cardiac output. Under the influence of a ramp input we see from Figure 1.18 that cardiac_ output rises and then begins to saturate as load increases. We also know that as load increases ve nous oxygen concentration falls as the exercising muscles extract more and more oxygen from the blood flowing through them. Thus, as the cardiac output (CO above) saturates, the venous oxygen concentration (Cvaz) falls. These comple ment each other in the above equation to produce the nearly linear output until at a high input level, the effect breaks down and the limiting effects of both CO and Cva2 appear and VOz plateaus. This complementary effect does not necessarily hold with other inputs and shows how important it is to use a variety of inputs in testing the system. In fact, what we see here is that the input "prepares" the venous pool in such a way that its decreasing oxygen concentration complements the increasing cardiac output to produce a nearly linear system output. Consider instead inputs that prepare the venous pool in different ways as shown in the top plate of Figure 2.4. 56
PAGE 64
150 (j) '"' '"' ""' "' 100 "' ""' 0 ,_, 50 40 40 u "" en ''"" 30 : '"' ., < 20 '"' o. ;:> "' 0 10 40 __________________________________ jc._ _______________ 1 60 80 100 120 TIME(SECS) 140 160 180 200 //___ j ___ / __ ,. 60 80 100 120 TIME(SECS) _j 140 160 180 200 Figure 2.4: Step up and step down inputs (top plate) and the resulting simulation outputs (bottom plates). The solid line is a step down and prepares the venous pool to be oxygen depleted. The dashed line is a step up and prepares the venous pool to be oxygen rich. The resulting outputs are shown in the bottom plate of Figure 2.4. If the system were truly linear the sum of the outputs would be the output of the sums (this is the definition of linearity). The top plate of Figure 2.5 shows the output resulting from the sum of the inputs and the bottom plate shows the sum of the outputs. Clearly they are different and the system is not linear. There are several reasons why the human 1>02 vs load system is not linear. 57
PAGE 65
48 u "' 47 UJ '.1 ::;: 46 t: ;:> 0.. 45 !:; 0 44 40 60 80 100 120 140 160 180 200 TIME(SECS) 48 u "' 47 UJ ',.., ::;: 46 t: ;:> 0.. 45 ... ;:> 0 44 40 60 80 100 120 140 160 180 200 TIME(SECS) Figure 2.5: The output of the simulation under the summed step inputs (top plate) and the sum of the outputs under each step input separately (bottom plate). 58
PAGE 66
100 (/) t: s 50 0 50 100 150 TIME(SECS) Figure 2.6: Sample step input. One is the complicated dynamics of the venous pool as discussed above. Another is that the time constants in the true human system are not exactly the same for the "off" dynamics as for the "on" dynamics and a third is that we allow the possibility of pushing the subject into the plateau portion of the curve in Figure 2.2 and this clearly produces nonlinear results. Further experimental results showing that the human system is nonlinear and a discussion can be found in [23]. Step Inputs The next input form we consider is the step input. A sample input is given m Figure 2.6 with the corresponding output in Figure 2.7. By averaging the output over many repetitions of real experimental data and by fitting curves to the data several investigators (for example, [23]) find that the human subject has the biphasic output of Figure 2.7. (Note in Figure 2.7 that for about fifteen 59
PAGE 67
30 (/) 20 E< z ::0 10 0 0 50 100 150 TIME(SECS) Figure 2. 7: Simulation under the step input. seconds after the step in the input that the shape of the output is different than the later portion giving the output a biphasic appearence.) Whipp et al. [46] speculate that the early phase consists of a rapid increase in cardiac output which is emptying a venous pool prepared by the low load portion of the step input to be oxygen rich. The second phase occurs as the venous pool becomes more and more oxygen depleted as a result of receiving oxygen depleted blood form the muscles responding to the high load portion of the step. Under a step input our simulation becomes equivalent to the model proposed by Barstow and Mole [5]. They show that their model is in agreement with experimental step data and so it is not surprising that our model also matches step data well. Figure 2.8 shows the simulation output with the output of an experimental subject superimposed. They agree well. In fact, we use this step data to determine the time constants of the linear portions of the underlying Wiener subsystems of Chapter One. Again, there are three degrees of freedomone for the linear subsystem in each of the 60
PAGE 68
20 __ / 10 0 50 100 150 TIME(SECS) Figure 2.8: Output of the simulation (solid) with output of a human subject (dashed) superimposed. Output units are ml. oxygen per second. three Wiener systems. We find that we can fit the data relatively well by choosing the three time constants to be the same and so we only use one degree of freedom in this fitting. This leaves two to model further subtleties in the data. To see how the choice of these time constants can effect the fine structure of the model consider the relation between the time constant for the Wiener subsystem that determines muscle blood flow (rqm) and the time constant for the Wiener subsystem that determines muscle oxygen consumption ( Tvm) If Tqm < Tvm an increase in load causes muscle blood flow to increase slower than muscle oxygen consumption and the effluent from the muscle falls in oxygen concentration as one would expect. However, if Tqm > Tvm, flow would increase faster than oxygen consumption and we would have the paradoxical result that the oxygen concentration of blood flowing out of the muscle would actually rise with increasing load. By varying the relative magnitudes of these two time constants 61
PAGE 69
150 tzl t: 100 z ;:> 50 0 20 40 60 80 100 120 140 160 180 200 TIME(SECS) 60 ... 40 /// / (/] / E< z ;:> 20 v 0 0 20 40 60 80 100 120 140 160 180 200 TIME(SECS) Figure 2.9: Some responses that can be obtained by fixing Tqm and varying Tvm. we can adjust the fine structure of the response to fit the data for a given subject. Figure 2.9 shows some of the variety of responses we can produce. Sine Wave Input The next input we consider is a sine wave. Figure 2.10 shows a .01 Hz. sine wave in the time domain in the top plate and in the frequency domain in the bottom plate. When this signal is sampled and held we get the more complicated signal shown in the time domain in Figure 2.11. In this example we will sample at ten second intervals. It is known [25] that sampling in the time domain has the 62
PAGE 70
200 150 (IJ ,.. ::;; "" 100 0 < 0 50 ,_, 0 0 50 100 150 200 250 300 TIME(SECS) 1 (IJ ,.. z :::> "' 0.5 "' "" 0 "" 0 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 FREQUENCY(HZ) Figure 2.10: Sine wave in the time domain (top) and frequency domain (bottom). Frequency spikes are at 0.01 Hz. and _Q.01 Hz. 150 (IJ ,.. 100 z :::> 50 50 100 150 TIME(SECS) 200 250 Figure 2.11: Sampled and held sinusoid in the time domain. 63 300
PAGE 71
PSD I I I FREOUENCY(Hz) Figure 2.12: Spectrum of the sinusoid after sampling. effect of replication of the spectrum in the frequency domain shifted successively by C..w = 21rjT, where Tis the time between samples.1 This means that the sinusoid will have the frequency domain representation shown in Figure 2.12. This sampled signal is then held constant for ten seconds until the next sample. This function is called a zero order hold (ZOH) (15] as opposed to a first order hold where two successive samples are linearly interpolated or a second order hold where successive sample are fit by a quadratic or other higher order holds. In the time domain this function can be represented as in Figure 2.13 which shows the impulse response of the ZOH. The transfer function of this system is the Fourier transform of the inipulse response and so is just the Fourier transform of a shifted 1This argument is often presented in conjunction with analyses of aliasing, where it is pointed out that if the samples in the time domain are too far apart (T too large) then Aw will be too small and the shifted copies of the spectrum will overlap resulting in aliasing. 64
PAGE 72
' INPUT TO ZOH (IMPULSE) .1 ZOHIITIME OUTPUT OF ZOH (STEP) Figure 2.13: Effect of ZOH on an impulse at t = 0. rectangular window: HzoH = exp( jwT /2)Tsinc(wT /2), TIME where sinc(x) = in;"l, if x f 0 and sinc(x) = 1 if x = 0. In terms of the power transfer function, we have: This is plotted in Figure 2.14. We can now cascade the sampling action with the hold action on the sinusoid by multiplication of the power transfer functions in the frequency domain to determine the output of the ZOH. Figure 2.15 shows this calculation by superimposing the graphs. From the point of view of testing the system the important point is that we are not presenting a sinusoid to The Bicycle subsystem and, thus, to The Subject subsystem, but a more complicated signal consisting of several frequency components. We would prefer to input a 65
PAGE 73
1 "'t 0.8 0.7 "' 0.6 t: z ;::> 0: 0.5 "' "" 0 0.4 0.. 0.3[ 0.2 0.1 0 5 4 3 2 1 0 1 2 3 5 FREQUENCY Figure 2.14: Power transfer function for the ZOH in normalized frequency units of (wT/2). 66
PAGE 74
Ul t: z ;:> "' "' "' 0 .,. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 L 0.2 0.1 0 5 / "'4 3 2 f\ I / \ / 1 0 1 2 3 4 5 FREQUENCY Figure 2.15: PSD of the sampled sinusoid superimposed on the power transfer function of the ZOH. Frequency is in units of wT /2. 67
PAGE 75
I _I 150 "' it: 100 z ;:> 50 50 100 150 TIME(SECS) 200 250 300 Figure 2.16: Sampled and held sinusoid after filtering through low pass filter with cutoff frequency of 0.1 Hz. sinusoid into the human system and the above suggests that we might do this by filtering. Recall that we still have the antialiasing filter to impose before passing the signal to the bike. If we were using this filter only to prevent aliasing the cutoff would be near the Nyquist frequency. If we are resampling (recall Chapter One) at, say, one sample every 0.5 seconds, the Nyquist frequency would be fN = 1Hz. However, we can certainly set it lower. Figure 2.16 is a plot of the output of the antialiasing filter if the cutoff is set at 0.1 Hz. If we set the cutoff frequency even lower we can finally remove all the frequency components except the fundamental as shown in Figure 2.17. Figure 2.18 shows both the original signal input to the computer and the resulting filtered output that is then input to the bike. We see that they are very similar except for a phase shift. Also note that the bicycle itself is a low pass filter. Recall from Chapter One that we model it as a linear system. Figure 2.19 shows the power transfer function for a 68
PAGE 76
200 150 l/l !::: 100 z p 50 0 0 50 100 150 200 250 300 TIME(SECS) Figure 2.17: Sampled and held sinusoid after filtering through low pass filter with cutoff frequency of 0.05 Hz. 200 .. .. / ... \ // \ / / 150 \ \ I I \ \\ / Ul I \ \ < / z 100 \ \\ \ p \ \ 50 \ // \, \ \ / '\ \ // '\ ,/ . / 0 .. 0 50 100 150 200 250 300 TIME(SECS) Figure 2.18: Input to the computer (solid) and final output to the bicycle (dashed). 69
PAGE 77
{/] < z ::> 0.5 0 200 150 100 50 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 FREQUENCY(HZ) Figure 2.19: Power transfer function for a typical bicycle. 50 100 150 TIME(SECS) 200 250 0.45 0.5 300 Figure 2.20: Output from passing the sampled sinusoid through the bicycle. typical bicycle and Figure 2.20 shows the result of passing the sinusoid through the bicycle alone without the antialiasing filter. It is clearly nearly sinusoidal. Thus, by proper filtering we can test the human system with a sinusoid in spite of the obligate sampling and holding that are part of the system. Figure 2.21 shows the result of testing a human with a 0.01 Hz. sinusoidal input with the output of the simulation superimposed. Clearly the two agree very well. Here we stress that in this case the simulation is not fitted to this data. The amplitude response of the simulation is determined by fitting the 70
PAGE 78
40 ' I ..... n J I I / tl 11 "' / If l 0 50 100 150 TIME(SECS) 200 250 300 Figure 2.21: Output of simulation and human subject under sinusoidal input. ramp response and the time constants of the simulation are determined by fitting the step response as discussed earlier. There is no fitting to the sine response yet they agree very well. The mean error (the mean of the quantity determined by subtracting the output predicted by the model from the observed value) for this study is 0.52 units. After we remove this bias the standard deviation of the residuals is 2.52 units. We will see later in this chapter that the noise has a standard deviation of 2.6 units about its mean and so this represents very good agreement. Figure 2.22 is the autocorrelation function of the residuals showing that they are nearly white. Figure 2.23 is the crosscorrelation function of the residuals with the input showing little cross correlation. All of these factors show that the simulation is doing a very good job of modeling the subject for the sinusoidal input data. 71
PAGE 79
1 0.5 t\ Ul ,... z p 0 0.5 0 2 4 6 8 10 12 14 16 18 LAG(SECS) Figure 2.22: Autocorrelation of the residuals for the sinusoidal input case. 1 0.5 Ul t: 0 z p 0.5 1 0 5 10 15 20 25 LAG{SECS) 30 35 40 45 20 50 Figure 2.23: Crosscorrelation of the residuals and the input of the sinusoidal input case. 72
PAGE 80
Constant Switching Pace Inputs In this section we consider the output of the simulation and the human under the input of a constant switching pace (CSP) signal. We will discuss this signal class more fully later and use them in the identification technique that will be developed. The signal is generated by selecting a hold time or viewed another way the time between switches from one level to another. During each hold time we choose the input level by a random draw from a uniform distribution. Figure 2.24 shows a CSP input and the resulting output of both the simulation and the human system. Again the outputs agree very well and as in the sinusoidal input case the simulation is not fitted to this specific signal. In this experiment the mean error was .43 and the standard deviation of the residuals was 2.69. The Model with Noise In this section we will introduce a simple n01se model into the simulation. Figure 2.25 shows a time domain plot of output from a human subject about its mean in the top plate and a frequency domain plot in the bottom plate under a constant input. The standard deviation of this signal is 2.6 units. Figure 2.26 shows a histogram of the amplitude distribution of this time series in the top plate and its autocorrelation function in the bottom plate. The histogram appears normally distributed. To verify this we plot in Figure 2.27 the residuals against 73
PAGE 81
30 100 (\ 1 \1\,1\t f 25 150 200 TIME(SECS) 200 TIME(SECS) 250 300 Figure 2.24: Constant switching pace input (top plate) and the corresponding output of the human system (dashed line, bottom plate) and the simulation (solid line, bottom plate). 74
PAGE 82
10 5 5 10 0 100 200 300 400 500 600 TIME(SECS) 1 0.5 "'
PAGE 83
150 <>: 100 "' :::> z 50 f0 8 6 4 .2 0 BIN 112 r 4 6 8 {J) t: z :::> "' 0 u 0 !:; .., 0 5 10 15 20 25 TIME LAG(SECS) Figure 2.26: Histogram of the noise amplitudes (top) and autocorrelation function of the noise (bottom). 76
PAGE 84
4 2 "' 0 0 u (I) z 2 4 6 4 2 0 2 4 6 VALUE Figure 2.27: Plot of residuals against the normal scores for noise data. the normal scores calculated by the MIN IT AB 2 routine NSCORE [33]. The plot is nearly linear as one would expect if the data were indeed normal. When we calculate the correlation coefficient of these data with a straight line we get .988 which allows us to accept the hypothesis that the data comes from a normal distribution by the ShapiroWilks test [33]. From the autocorrelation function of Figure 2.26 we see that the noise is nearly white. Thus, at least to a first approximation the output noise is nearly white and nearly Gaussian. The near Gaussianity can be explained by assuming that the final noise results as the sum of many small, independent sources and then applying the central limit theorem to see that the final noise is near normal. In order to understand the near whiteness of the noise we observe that the cascaded subsystems are smoothing and so noise effects tend to be removed as data flows through the system. Each Wiener subsystem consists of a linear system which acts as a low pass filter. 2 Registered trademark of Mini tab, Inc. 77
PAGE 85
(I] f" ';;; i!= "' .., 0 .., u "" (I] '.., ::;: ;::;::> 0.. f" ::> 0 200 150 100 50 ;u 0 0 60 40 20 0 0 "'""\ m; 50 100 150 50 100 150 II ll JJ rM ru \'" 200 250 300 350 400 450 TIME(SECS) 200 250 300 350 400 450 TIME(SECS) Figure 2.28: Noisefree input (top) and corresponding output (bottom). 500 500 Similarly, the venous pool acts as a smoothing system. It averages the inputs with the previous state of the pool to produce the current oxygen concentration in the pool. All of these effects tend to remove the noise effects from previous stages. This is illustrated in Figure 2.28 which shows the result (bottom plate) of inputing a noisefree signal (top plate) into the simulation. Figure 2.29 shows the same input (top plate) corrupted with Gaussian white input noise with power level at ten percent of the original input. The output, shown in the bottom plate, is essentially unchanged from the output in the noisefree case. Thus, we will in this thesis use additive, white Gaussian output noise as our first approximation 78
PAGE 86
300 u "' (IJ 200 'w M ..., A :>: ,:;:.:> 100 0.. E< :.:> 0 0 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) 60 u "' (IJ 40 '..., :>: ,:;:.:> 20 0.. E< :.:> 0 0 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) Figure 2.29: Noisy input (top) and corresponding output (bottom). 79
PAGE 87
to a noise model in the simulation. As more data on the exact values of the noise in each subsystem become available we will incorporate it. Summary In this chapter we have studied the output of the simulation and the output of human subjects under several standard inputs and introduced a noise model into the simulation. The results show very good agreement between the simulation and experimental data. An exhaustive validation would require fitting the simulation to many subjects under a variety of experimental conditions. We simply do not have the resources available to do this and will leave this for future work. However, we feel that the physiology we have used to develop the simulation is correct and that the examples presented here show that the simulation is an excellent first order approximation to the human system. We also feel that its incorporation of simple and complex nonlinear behaviors makes it very useful as a testbed for the development of system identification algorithms. In the next chapter we will develop such an algorithm and in Chapter Four we will apply it to the simulation. 80
PAGE 88
CHAPTER 3 DEVELOPMENT OF THE IDENTIFICATION METHOD Introduction In this chapter we will begin the process of constructing a method of system identification. We here again state our objective of determining, by using only input and output data, a global identification in terms of polynomial functions of lagged inputs and outputs. If we were dealing with human subjects the objective would be to design an experimental protocol that would be run on a given subject to produce a mathematical model of the subject. Such a mathematical model would be ideal if it produced the same output as the subject for every possible input. Practically, such an ideal model can not be constructed and we instead construct models that approximately match the subject output for some limited class of input signals. In this thesis we will instead be identifying the computer model of the human subject developed in Chapter One. Although the identification of this relatively complicated model is interesting in its own right, our idea is to use it to determine and solve problems associated with identifying the human system. The problems that arise are parallel in both systems; however, the simulated system allows us much more control. For example, in the computer simulation we know exactly the structure of the system, its parameters 81
PAGE 89
and its noise sources. We also have access to internal variables not available in human experiments and we have the ability to alter all of these quickly and easily without putting human subjects at risk or exposing them to long, difficult experimental sessions. Because the techniques we develop relate to both the simulated and real systems, we will freely mix notation in the following and say, for example, "the input to the bicycle is passed to the subject" when, in fact, we mean "the input to the simulated bicycle is passed to the simulated subject". If there are situations where this association of human based experiment and simulation based experiment break down we will point them out. In developing the identification there are several questions that must be addressed. According to Ljung [29] we must determine three entitiesthe data record, the set of candidate models and the rule for assessing the models. Ineluded in the data record are such questions as the form of the input, what outputs to measure, when to sample both the input and the output and other questions related to data collection. The candidate models are determined by how we represent the system and the assessment rule is the criterion used to determine which model in the candidate set "best" represents the system. We must also establish an identification method that selects from the candidate models the ones that optimize the assessment rule. In the first part of this chapter we will discuss each of these aspects in turn. In the second part of the chapter we will apply them to the identification of a simple system in preparation for applying them to the full simulated system in Chapter Four. 82
PAGE 90
The Data Record The data record we consider is an inputoutput pair where the input is the test input the investigator enters into the computer and the output is an interpolated and smoothed version of the breathtobreath oxygen consumption data as discussed in Chapter One. Regarding the input data, there are two constraints imposed by the experimental system. One is that the input level is limited so that the work load presented to the subject is between zero watts and about 250 watts. While it would be possible to design a driven bicycle and a coupling mechanism so that energy transfer would be from bike to subject, this is clearly a different system than we are trying to test and so the lower limit is unloaded cyclingzero watts. At the other extreme, there is a limit above which the subject cannot maintain constant pedaling speed on the bike and so we must limit the upper range of the input to about 250 watts. The second constraint is that the duration of the experiment must be no longer than about thirty minutes. This is because subjects fatigue and if the experiment is too long physiologic parameters change and the system is no longer stationary. Both of the above constraints are subject dependentan elderly patient with emphysema will have different limits than a young endurance athlete. The ranges considered here are for healthy sub jects and we will discuss the effects of limiting both load range and experiment duration. Another limitation on the input anses from the experimental setup as dis83
PAGE 91
Figure 3.1: The general modeling scheme. cussed in Chapter One. Although it would be possible to drive the bicycle with a continuous analog signal, the .typical experimental setup consists of sampling the continuous input at a constant rate and presenting this sample to the bike for a constant hold time equal to the interval between samples. This technique of sampiing followed by a zero order hold is an artifact of the experimental hardware, but is also a common technique in control theory [2]. Figure 3.1 shows the overall modeling scheme. Arbitrary continuous inputs (A) pass through the "Ideal system" to produce the corresponding continuous outputs (a). In reality, since the experimental system can pass only piecewise constant inputs to the bicycle, the actual experimental loop is (A) to (B) to (b) to (a). Thus, the arbitrary in84
PAGE 92
puts at (A) pass through the Sample and Hold to become the piecewise constant inputs at (B). These then are presented to the "True System". Inside the "True System" these inputs produce the breathtobreath data that are interpolated as discussed in Chapter One to produce the continuous outputs at (b). These outputs are equivalent to the output of the "Ideal System" in the sense that the distance between them in the chosen norm is less than some specified tolerance. It is this tolerance that will determine the upper bound on the hold time: If the hold time is too long, the signals at (B) will be too different from the signals at (A) and the resulting signals at (a) and (b) will differ by more than the required tolerance. Figure 3.2 shows a .01 Hz. sine wave and the resulting sample and hold signal for hold times of 5 seconds (top plate) and 20 seconds (bottom plate). It is clear that the 20 second hold signal is less like the original signal than the 5 second hold version. Therefore, we would expect that a 5 second hold in the sample and hold of Figure 3.1 would produce a signal at point (a) more like the result produced by the "Ideal System" than would a 20 second hold. The question of whether they are within the specified tolerance is one we will quantify later in this chapter when we consider a specific system. There are many different possibilities for determining the hold time structure in the identification experiment. One method is to let it be a random variable as discussed in Tulleken [43]. In this thesis we will require that the hold time be a constant. The determination of this constant is based on the general tradeoff 85
PAGE 93
1 J ' ' [fJ E< z :::> 0.5 : 1 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) 1 ' , r i i i i 0.5 l ...! ' I i ' [fJ \ l 1 t: i 0 ___ j I z \ I 1 :::> i 1 I 0.5 1 , , , , I 't ' i I I 1 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) Figure 3.2: A .01 Hz. sme wave and its sample and hold for hold times of 5 seconds (upper) and 20 seconds (lower). 86
PAGE 94
that if the hold time is too long, then signals at points (a) and (b) in Figure 3.1 do not agree within the tolerance discussed above. On the other hand, if the hold time is too short then the input will be deficient in low frequencies and will not stimulate important modes of the system and will not allow adequate excursion in the amplitude of the signals input to the nonlinearities to allow complete identification. The other factor we must consider is the amplitude during each hold time. If the input is a sine wave or other deterministic signal the amplitude during the hold time is the value of the input at the sample time. In designing test inputs we will take a slightly different point of view (recall the discussion in The Overall System section of Chapter One) and choose the amplitude value during a given hold time as a draw from a probability distribution. Several possibilities have been studied. Marmarelis [30] considers, among other possibilities, independent, Gaussian values. Barker [4] and many others consider values that alternate between two values according to a very specific algorithm which produces binary pseudorandom values. In this thesis we will follow Billings and Voon [8] and use independent, uniformly distributed values. In the analysis of a linear system or in the local analysis of nonlinear systems there is much to be said for pseudorandom sequences. By construction they have peaked autocorrelation functions (and therefore broad power spectral densities) that are very useful in identifi.cation [13]. Their shortcoming in the global identification of nonlinear systems is illustrated in Figure 3.3. Each curve represents the inputoutput function of a 87
PAGE 95
100 (/) E< "' :::> E< 50 :::> "" E< :::> 0 0 0 10 20 30 40 50 60 70 80 90 100 INPUT UNITS 100 [/) '"' z :::> '"' 50 :::> "" '"' :::> 0 0 0 10 20 30 40 50 60 70 80 90 100 INPUT UNITS Figure 3.3: Two nonlinear inputoutput functions. zeromemory nonlinear system. This means that for a given input on the xaxis the output is the corresponding y value. The two curves have the same y values at the x values of 25 and 75. If we were to choose these as the values for the two levels the systems would appear the same. Just as a line is determined by two points such bilevel methods work for linear systems, however, for nonlinear systems they do not work in general. Gaussian distributed values also have a similar shortcoming as illustrated in Figure 3.4 which shows a zero memory nonlinear transfer function with two Gaussian distributions superimposed. One way to view the identification process is that we approximate the nonlinearity with polynomial terms in such a way as to minimize errors at points along the xaxis 88
PAGE 96
5 4.5 4 3.5 {/) t: 3 z ;::> ,... 2.5 ;::> '" ,... 2 ;::> 0 1.5 1 0.5 0 1 2 3 4 5 6 7 8 9 INPUT UNITS Figure 3.4: Nonlinear inputoutput function with superimposed gaussian distributions. 89
PAGE 97
drawn from the distributions. Clearly the coefficients of the polynomials and, therefore, the model resulting from identification will be different for the two different sampling distributions. Since we hope to find a model valid across the entire range of amplitudes the uniform distribution is an obvious choice. In the nonlinear systems we consider, which are cascades of linear and zeromemory nonlinear systems, there are other ways to insure adequate excursions in amplitude at the input to the non linearity, but in this thesis we will use the uniform distribution and the amplitude values during the hold times will be uniformly distributed across the range of the input. Thus, the test inputs that we will use for creating the data record will be signals that switch values at multiples of a fundamental time called the hold time. During each hold time the value will be determined by a sample from a uniform distribution between the lowest possible input value and the highest possible input value. An example of such a signal is given in the top plate of Figure 3.5 for a 5 second hold time. Such signals are often called staircase signals because of their appearance, but we will call them constant switching pace (CSP) signals [30]. In the actual identification we will use methods that involve the crossand autocovariance properties of the inputs and outputs. One useful property is for the autocovariance function of the input to approach a delta function. A generalization of this is a concept which deals with the condition number and invertability of the matrix of lagged autocorrelation functions called persistent excitation (see Ljung [29]). 90
PAGE 98
Ul z ::> 0.5 0 0 50 100 150 200 250 300 TIME(SECS) 1 0.5 30 20 10 0 LAG 10 20 30 Figure 3.5: Constant switching pace signal with 5 second hold time (top plate) and its autocovariance function (bottom plate). It is easy to show [41] that the autocovariance function for the inputs we are using ( CSP signals) is NT 0 <:;_ T <:;_ N rr r(r) = N+r N <:;_ T <:;_ 0 N (3.1) 0, otherwise. where N is the hold time and that our input is persistently exciting of all orders. This autocovariance function is plotted in Figure 3.5 (bottom plate) for N=5 and consists of a triangular pulse. It is well known that as such functions become more localized in time their Fourier transforms spread in frequency. This will 91
PAGE 99
allow us to match the bandwidth of the input to the bandwidth of the system as we will discuss in detail when we consider the specific system. There are two further points we must consider. One is that the experimental conditions require that the identification be done in a limited time and so we are not really dealing with infinite sequences. This raises the question of how many samples are needed for the finite sequences to approximate the properties of the infinite sequences closely enough to be of use and the related question of determining if some finite realizations are better than others in the practical identification. The other point is that in linear system identification we need only know the properties of the inputs and outputs up to their second momentsie. means and covariances or corresponding spectra. In nonlinear systems we must, in general, know higher order moments or polyspectra for a complete analysis [35]. This raises the question of the higher order properties of the inputs we are using. Neither of the above points appear directly in the identification method we will develop and so we will not discuss them further in this thesis. The question of the best finite realization, however, is quite important practically and we will include it as future work. To complete the discussion of the data record consider the sampling of the output during the identification experiment. It this thesis we will sample the output synchronously with the input. There are many other possibilities as discussed in Crochiere and Rabiner [12]. In nonlinear systems there may be some advantage to sampling the output faster than the input because harmonics of 92
PAGE 100
the input frequencies can appear. We choose, however, to keep the processing bandwidth the same for both input and output and deal with this problem by widening the overall processing bandwidth. This will be discussed further in the setting of the specific systems. Candidate Models There are several possible candidate mathematical models we could use [1], [3]. There are also several ways to deal with the problem of identifying a naturally continuous system [37], [44]. However, in this thesis we will sample both the input and output processes and construct a discrete model between these sampled processes. Once this discrete model is identified, the way the overall model will work is by sampling the given continuous input, passing these samples through the discrete model to produce output samples and then these will be interpolated to produce the final continuous output. This procedure is shown by the lower loop in Figure 3.1. Arbitrary continuous inputs (A) pass through the "Ideal System" to produce the corresponding continuous outputs (a). The corresponding model loop will be (A) to (B) to (C) to (c) to (b) to (a). Although this process may appear complicated it is commonly used in both system identification and control applications. The principal advantage being that the discrete model can be determined using digital computing techniques. In this thesis we will restrict ourselves to two general types of models. One is the Volterra series as presented in Rugh [36] and the second is the Nonlinear 93
PAGE 101
AutoRegressive Moving Average with eXogenous inputs (NARMAX) model as presented in Leontaritis and Billings [27]. The Volterra model represents the output as a polynomial function of the previous inputs. The N ARM AX model represents the output as a nonlinear function of previous inputs and outputs. If we further approximate this function with a polynomial in its variables we get the polynomial N ARM AX model. The Volterra Series A sum of the form (3.2) is called a Volterra series and the hn are called Volterra kernels. The expansion was introduced by Volterra around 1900 as a generalization of power series to functionals and has been called a power series with memory [36]. If the series converges for a given set of time functions X = { x( t )}, then the series defines a functional on X. This functional can be thought of as mapping the "past" of xfrom time zero to time t into the current value of yy( t ). As such it is a representation of an inputoutput system with x(t) as input and y(t) as output. There are many more or less mathematical justifications for this method of representing systems discussed in the books by Schetzen and Rugh [39], [36]. We note here that both the Hammerstein and Wiener models discussed in Chapter One fit into this representation when the nonlinearity can be represented 94
PAGE 102
as a Taylor series. Recall that both of these models consisted of a linear system that contains the dynamics and a zeromemory nonlinearity. For example, let the nonlinearity be where z is the input to the nonlinear system and N(z) is the corresponding output. We can represent the linear part of the system as a convolution [7] L(t) = (oo h(tr)x(r)dr, lo where x is the input to the linear system and h is the impulse response. Consider first the Wiener model (Recall Figure 1.7) in which the input first goes into the linear system and the output of the linear system then becomes the input to the nonlinear system. Let the input be x(t) and so the output of the linear part of the system, after a change of variable and using causality and the fact that x( t) is zero for t < 0, is given by L(t) = f" h(r)x(tr)dr. This becomes the input to the nonlinearity and so we substitute it into the polynomial representation above to get the output of the nonlinear system, which is also the output of the overall system y(t) = N(L(t)) = a0 + a1L(t) + a2(L(t))2 + We then substitute the integral representation of L into this to get y(t) = ao +at j h( r)x(tr)dr + a2 j h(rt)x(trl)dr1 j h( r2)x(tr2)dr2 + ... 95
PAGE 103
Interchanging the order of integration we find that This is in the form of a Volterra series with kernels h1(r) = h(r), h2(rhr2 ) = h(r1)h(r2), and hn = h(r1 ) h(rn) This shows that the Wiener model has a corresponding Volterra representation. Now consider a Hammerstein model (Figure 1.7) and recall that here the zero memory nonlinearity comes first followed by the linear subsystem. Again let :v( t) be the input to the overall system and so the output of the nonlinear subsystem is obtained by putting this value into the above polynomial to get N(:v(t)) = ao + a1:v(t) + a2(:v(t))2 + This output of the nonlinear subsystem then becomes the input to the linear subsystem and so the overall output of the system is obtained by putting N(:v(t)) into the convolution representation of the linear subsystem. This gives y(t) =I h(r)N(:v(tr))dr. Using the polynomial representation of the nonlinearity yields y(t) =I h(r)[ao + a 1:v(tr) + a2(:v(tr))2 + ]dr. If we then distribute the integral we get, finally, that the overall system output IS y(t) =I a0h(r)dr +I a1h(r):v(tr)dr +I a2h(r)(:v(tr))2dr + 96
PAGE 104
Again, this has the form of a Volterra senes when we let the kernels be rz), and Since we will be building only discrete time models we now discretize the above input'output Volterra representation. The inputs, as discussed earlier, are piecewise constant because of the nature of the testing equipment, we can let the input value in the interval [iS, ( i + 1 )S). where S is the hold time, be x( i). We will further assume that the system has a finite memory of N S, where N is an integer. This means that each h,. "" 0 if any of its argumentsrqare greater than N S or, equivalently, that the system has a settling time approximately equal to N S. Omitting the S for notational simplicity and decomposing the integrals in equation 3.1 into regions of constant input we get that N N N y(n) = ho + L H1(i)x(ni) + L L H2(i,j)x(ni)x(nj) + (3.3) where etc. If we assume that the h's are smooth and make S sufficiently small that h" h2 vary little over [iS,(i + 1)5) it follows that H1(i)"" Sh1(i), H2(i,j)"" S2h2(i,j), Written in this way we see that equation 3.3 is a representation of the output of the system in terms of polynomials in the lagged input to the system. This is the Volterra series representation of the system. 97
PAGE 105
The NARMAX Model Goodwin and Payne [18] develop several parameter estimation algorithms for linear systems based on representing the system as a linear difference equation: n m y(t) =I:a;y(ti) + L b;x(ti), i=l i=l where y(t) is the output at timet, y(ti) is the output at a time i*(sampling time) prior tot and x(ti) is the input i sampling times prior tot. The nonlinear analogue of this is y(t) = F[y( t1 ), y( t 2), y(tn ), x( t1 ), x( tm )], where y and x are as above and F is a nonlinear function. Leontartis and Billings [27] develop this model under the name N ARM AX for nonlinear autoregressive moving average with exogenous inputs and show that the Hammerstein and Wiener models are special cases of this representation. It is also clear that the Volterra model of the previous section is a special case if we suppress the arguments that involve lagged outputs. Part of their justification for such a model is that since it includes output, which depends recursively on previous inputs, the NARMAX method may yield a more parsimonious model (that is, one with fewer parameters) than the Volterra model. Continuing along the lines of the linear ARMAX model [19] we add a noise model to this system. Following Billings and Voon [8], we assume that the noise is additive at the system output and so the observed output will be z( t) = y( t) + e( t ). We solve this for y to get that 98
PAGE 106
y(t) = z(t)e(t) and put this in the original equation to get the model z(t) = F[z(t 1), z(tn), x(t 1), x(tm), e(t 1), e(t I)]+ e(t) where we note that the e can not be factored out into an additiveonly term as it can be in the linear case, but is included in the structure of F. Further, if P can be approximated by a polynomial in its arguments we get the final polynomial regression form that we will use as our second possible model form. As an example of how such a NARMAX model might arise suppose the underlying system is a Wiener system consisting of a first order linear system corresponding to the differential equation ry + y = x and the subsequent zero memory nonlinearity (recall Figure 1. 7) is a squarer represented by z( t) = [y( t )j2. If we discretize the linear subsystem in small time steps of magnitude h we get y(t+h)y(t) () () T h +yt ""'xt. We can solve this for y( t + h) in terms of previous inputs and outputs to get h h y(t +h)= (1)y(t) + x(t). T T This is then the input to the zero memory nonlinearity (the squarer) and so the final system output is h. h h h z(t +h)= (1j2(y(t))2 + ()2(x(t))2 + 2(1)()y(t)x(t). 7 T 7 T This is a N ARM AX model for this Wiener system. 99
PAGE 107
The Assessment Rule In this section we determine how to evaluate the models and how we will determine the parameters in the model. The simple linear regression model is written where y(t) is the system output at timet, if is a parameter vector (we will use the vector notation to mean column vectors) and J;( t) is a vector of regressors [17]. An experimental procedure produces observations of y and ;j at various times and the identification problem is to find a if so that the model in some sense represents the observations. One tries to design the experiment so that the observations are rich enough that the model will represent the output y for other observations also. For multiple observations the above is simply a system of equations and aswe get more and more observations the system becomes overdetermined. One way to solve the above set of equations is to define the output error to be e(t) = y(t);f(t)if. Then we construct the vector of e's for all observations and choose if to minimize e in some vector norm. If we choose the [2 norm the resulting estimation of if is called the least squares estimate. If we index the observations we can construct the matrix equation 100
PAGE 108
where y(l) y= y(N) and ,V(ll P= jT(N) and e is the vector of errors. The vector of errors in this setting can be minimized analytically [45] and the solution is This is called the least squares estimate of if. Under the assumption that each component of e is independent with zero mean and uncorrelated with the regressors, it is known that B is an unbiased estimate of B. Furthermore, among estimators of the form (linear estimators), the one above has the smallest covariance matrix in the sense that for any Q, Q(PTP)lpT is nonnegative definite. Thus, B is called the best linear unbiased estimate of B. In addition, if each of the e( t) 's is Gaussian, then B corresponds to the maximum likelihood estimate of B. Now note that if we let ;fY(t) = (x(tl),x(t2), .. ,x(tn), 101
PAGE 109
a?(t1),x(t 1)x(t2), ,x(tn)2 x3(t1),), then the form of the regression equation agrees with the form of the Volterra series and the components of if are estimates of the Volterra kernels. Also, if we include in ;f(t) polynomials in the lagged output, the form of the regression equation agrees with the polynomial NARMAX model. In this thesis we will use these forms to construct models of our system and the assessment method will be to minimize e or some function of it. We will begin with simple least squares with the regressors being the lagged inputs and/ or lagged outputs as above. This will lead to difficulties as noise is introduced into the model. It is not surprising that noise would lead to increased variance in the estimates, however, noise can also lead to bias. To see why this is true, we follow Soderstrom and Stoica [41] and recall that Suppose that the data are generated by the noisy process y(t) = ;fif + n(t), where if is the "true" parameter and n IS nmse. In the matrix notation this becomes i/=il>O+N. 102
PAGE 110
Thus we can write the difference between the "true" and estimated parameter vector as iJe = (TtlTy () (Ttl[Ty ( T )Bj = (TflT[y Bj = ( T f 1 T [NJ. Thus, in general, if and N are correlated the estimate of () will be biased. One advantage of using a simulation model is that we can determine the effect of different noise models and signal to noise ratios on the quality of the estimation. There have been methods developed to estimate () and we will discuss them and determine their effect in Chapter Four. An Example Before we begin the analysis on the full simulation, we will illustrate the preceeding discussion with a simpler nonlinear system. Let S be a Wiener system (recall Figure 1. 7) in which the linear subsystem L is a first order linear system with r = 10 sees. Thus, the differential equation governing this linear system is lOy+ y = X. Also let the nonlinear subsystem be a "clipper" with the inputoutput function shown in Figure 3.6. Suppose that the range of inputs to the overall system is 103
PAGE 111
1 0.9 0.8 0.7 Ul .... 0.6 z :::> .... 0.5 :::> 0.. .... 0.4 :::> 0 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 INPUT UNITS Figure 3.6: Inputoutput function for the clipper. 104
PAGE 112
1 Ul 0.5 ::0 0 0 \ \ \ 100 f\ I\ \i \I \ 200 300 400 I I 1\ I I 1\ (\ \; \I \i \I \/ \ 500 600 700 600 900 1000 TIME(SECS) 12 z 0.5 ::0 100 200 300 400 500 600 700 600 900 1000 TIME(SECS) Figure 3.7: Sine wave input (top) and output (bottom) of the nonlinear system. between zero and one. This system is much like the Wiener models we use in constructing the first stages of the model in Chapter One; the major difference is that the clipper is a more dramatic nonlinearity than the MichaelisMenton type curves used there. If we input a sine wave of frequency 0.01 Hz. as in the top part of Figure 3. 7 into the system we get the output as a "clipped" sine wave of the same frequency as shown in the bottom part of Figure 3. 7. This is expected because the linear system produces a sine wave of the same frequency as the input [7]. This then becomes the input to the nonlinear subsystem which then simply "clips" off the tops to produce the output. Now we input a constant switching 105
PAGE 113
pace (CSP) signal into the system. This is the type of input discussed earlier where the amplitude remains constant between switching times every S seconds and the amplitude during each interval is drawn from a uniform distribution between zero and one. Figure 3.8 (top plate) shows such an input. In this case the hold time is only one second and so does not stand out in the figure. Under this input the linear system produces its standard step response as it traverses between levels to produce the output in the center plate. The zero memory clipper then "clips" the tops off to produce the final system output. To generate the input we will use for the identification experiment we generate samples from a uniform distribution between zero and one as discussed in The Data Record section. These will be the input values between the switch times. We then match the bandwidth of the input to the system so that input power is not wasted stimulating modes of the system that do not respond. This is done as discussed previously by sampling and holding the input. Thus, the first step is to determine which frequencies are passed by the system. This can be done analytically as follows. First note that in the frequency domain, power is attenuated by the linear subsystem according to 1 P(f) = (21rjr)2 + 1' where f is the frequency in Hertz and r is the time constant of the linear systern (10 seconds in this case). This is derived from the equation of the linear 106
PAGE 114
(/J ,.. z 0.5 p 200 300 TIME(SECS) 400 500 600 (/J ,.. z p (/J !::: z p 0.5 ',!' 0 0 100 0.5 200 300 TIME(SECS) 400 500 600 OL______ _L ________ L_ ______ _L ________ 0 100 200 300 TIME(SECS) 400 500 600 Figure 3.8: Constant switching pace input (top), output of the linear subsystem (middle) and output of the nonlinear system (bottom) for the "clipper" system. 107
PAGE 115
subsystem: riJ + y = x. In the Laplace domain, this equation becomes rsY + Y =X, where s is the Laplace domain variable, capitals represent Laplace transforms and we have chosen y(O) = 0. Thus, Y(s) 1 X(s) rs+1 Converting to the angular frequency domain (s = jw) we get Y(w) 1 X(w) = rjw + 1 The corresponding power transfer function is Converting this from radians per second to Hertz by letting w = 27r f we get 1 P(f) = (27rfr)2 + 1' where f is the frequency in Hertz. A plot of this power transfer function is shown in Figure 3.9. From this formula or the plot we see that the power decreases to of its maximum at f = .0159 Hz. for this linear system with r = 10 sees. In order to determine the effect of the nonlinear clipper in the frequency domain we input a single sine wave into the system. We choose a frequency of .01 Hz. This 108
PAGE 116
1 0.9 0.8 0.7 z 0.6 0 E::: .., ;:o z 0.5 "' 0.4 .., 0.3 0.2 0.1 0 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 FREQ(HZ) Figure 3.9: Plot of power attenuation by the linear subsystem vs frequency. 109
PAGE 117
1 '\ (\ (\ f\ f\ 1\ 1\ (\ 1\ 1\ [/) t:: 0.5 z ::> 0 \i \/ \; \I \; \) \/ \/ \J \J 0 100 200 300 400 500 600 700 BOO 900 1000 TIME(SECS) 1 [/) t:: z 0.5 ::> 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 FREQ(HZ) Figure 3.10: A .01 Hz sine wave in the time domain (top) and frequency domain (bottom). 110
PAGE 118
10 8 6 4 2 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 FREQ(HZ) Figure 3.11: Frequency domain output from sine wave input into the clipper. signal is shown in both the time domain (top) and frequency domain (bottom) in Figure 3.10. The corresponding frequency domain output is shown in Figure 3.11. We know that the linear subsystem will output a sine wave at the same frequency (0.01 Hz.) as the input. This then becomes the input to the nonlinearity and, thus, the clipper, as we can see from Figure 3.11 produces a significant response at the fundamental frequency and at the first and third harmonics corresponding to frequencies of J, 2f and 4f respectively. As we have argued (Appendix A), this allows us to write the overall system output as z(t) = A.1 exp(j2rrft) + A2 exp(j4rrft) + A1 exp(8rrft), 111
PAGE 119
I I which corresponds to expanding the nonlinearity as a polynomial The coefficients A; decrease rapidly (as can be seen by noting the magnitude of the spectral components in Figure 3.11) and so we feel safe in saying that the overaJl system output will contain almost no power above 4*(bandwidth) = .0636 Hz. Therefore, if we allow for frequencies up to .1 Hz. we can be quite certain that we are not missing frequency domain information. The Nyquist frequency corresponding to .1 Hz. is .2 Hz., ie. one sample every 5 seconds. We know the value of oversampling but when we count the dual effects of the limited bandwidth of the linear subsystem followed by the rapid decrease of the coeffients of the nonlinearity, we feel that we are not missing significant high frequency information either in the input or the output if we sample every 5 seconds. In the earlier section The Data Record we showed that the autocovariance function of the CSP inputs we are using was a triangular pulse. We also noted that narrowing this pulse in the time domain spreads its representation in the frequency domain and conversely. We also noted that as we shorten the hold time of the CSP signal we narrow the triangular pulse of its autocovariance function and conversely. Thus, we have a way to match the spectrum of the input with the bandwidth of the system we are trying to identify. To determine the spectral content of a CSP input with hold time of N seconds we Fourier transform its 112
PAGE 120
1 10 Ul Ul '"' 0.5 t: 5 z z ::0 ::0 0 0 0 50 0 0.5 LAG FREQ 1 Ul Ul '"' 0.5 t: 5 z z ::0 ::0 0 0 0 50 0 0.5 LAG FREQ Figure 3.12: Spectrum (right plates) corresponding to the triangular autocorrelation (left plates) for N =2 (upper set )and N = 10 (lower set). autocorrelation function ( eqn. 3.1) to get FT(w) =_!._sin( l. N sin(}) Figure 3.12 shows the autocorrelation function for N = 2 and N = 10 in the upper and lower left hand plates and the corresponding Fourier transforms in the right hand plate. Note that as N increases, FT narrows. Clearly we can affect the bandwidth of the input by changing the hold times. To see this effect more clearly in the setting of our system, Figure 3.13 shows a time domain input consisting of a constant switching pace signal with a one second hold in the upper figure. 113
PAGE 121
The middle figure is the output of the linear subsystem and the lower figure is the final system output. The corresponding power spectral densities are shown in Figure 3.14. We see two phenomena in these figures. One is that we are needlessly wasting power in the high frequencies that are input to the linear subsystem, but are not to the nonlinear subsystem. The other is that although the input to the linear subsystem covers the amplitude range quite well, the corresponding output of the linear subsystem (which is the input to the nonlinear subsystem) is quite restricted. The series of Figures 3.15 to 3.18 shows similar plots for holding times of 2 seconds, 5 seconds, 10 seconds and 20 seconds. In these figures the upper left plot is the input verses time and the lower left is the output of the linear subsystem verses time. (We plot this instead of the output of the clipper because we want to see the excursion of the input to the nonlinear system. The output of the total system is just this signal with the tops cut off.) The upper right figure is the power spectral density (PSD) of the input calculated by Welch's method [26] and the lower right figure is the PSD of the total system output. We note that the 5 second signal satisfies the requirement that the input band extends beyond the output band and so there is no loss of frequency information and yet there is not a lot of power expended on frequencies to which the system will not respond. Figure 3.19 shows in its top plate a histogram of the amplitude samples for the input to the system for the one second hold experiment. It is quite flat and the entire range is well represented. The successive plates in Figure 3.19 show 114
PAGE 122
UJ ,.. 0.5 z "' 0 0 100 200 300 400 500 600 TIME(SECS) 1 til ,.. z "' 0 0 100 zoo 300 400 500 600 TIME(SECS) 1 til ,.. z 0.5 "' 0 100 200 300 TIME(SECS) 400 500 Figure 3.13: Time domain signals for 1 second hold input. 115 600
PAGE 123
UJ E< z 0.5 ::> 0.5 0.05 O A 0 0.05 z 0.5 ::> 0 0.05 0.1 0.15 0.1 0.15 0.1 0.15 0.2 0.2 FREQ(HZ) 0.25 FREQ(HZ) 0.25 FREQ(HZ) 0.3 0.3 0.35 0.4 0.45 0.35 0.4 0.45 Figure 3.14: Frequency domain signals for 1 second hold input. 116 0.5 0.5
PAGE 124
"' ,... z 0.5 p 500 TIME(SECS) 1,, 0 500 TIME(SECS) 1,, 0.5 0.5 FREQ(HZ) 1,, f!.; z 0.5 p 0.5 FREQ(HZ) Figure 3.15: Time domain input (top left) and output (bottom left) and corresponding frequency domain input (top left) and output (bottom right) for a 2 second hold time. 117
PAGE 125
Ul Ul t: 0.5 t: 0.5 z z ::::> ::::> 0 0 0 500 0 0.5 TIME(SECS) FREQ(HZ) 1 1 Ul Ul t: 0.5 t: 0.5 z z ::::> ::::> 0 0 500 0 0.5 TIME(SECS) FREQ(HZ) Figure 3.16: Same as previous figure except for 5 second hold time. 118
PAGE 126
1 "' 0.5 0 0 J llrJ IJ 500 TIME(SECS) 1., z 0.5 \:l 1.. N 0.51>, \ V\. \ 0 0.5 FREQ(HZ) 1,0.5R 0 500 0 0.5 TIME(SECS) FREQ(HZ) Figure 3.17: Same as previous figure except for 10 second hold time. 119
PAGE 127
1 1 !:1 I (/] (/] t: 0.5 f< 0.5 z z I ;::> ;::> i \ \ \ 0 \A. 500 0 0.5 TIME(SECS) FREQ(HZ) 1 1 (fJ "' f< t: 0.5 z z ;::> ;::> 0 0 \ 0 500 0 0.5 TIME(SECS) FREQ(HZ) Figure 3.18: Same as previous figure except for 20 second hold time. 120
PAGE 128
the histograms of the amplitudes of the output of the linear subsystem (which is the input to the nonlinear subsystem) for hold times of one second, five seconds and ten seconds respectively. Here we see that as the hold times increase the amplitude excursion of the output of the linear system increases. This can also be seen in the previous series of figuresduring each hold time the first order linear subsystem produces its step response and the longer the hold time the larger the excursion. This is the second part of the trade off we mentioned earlier in this chapter. When the hold time is short the output of the linear subsystem is a poor test input for the global identification of the nonlinear subsystem because it is clustered around the output mean. The compromise we reach is the 5 second hold signal, which is a reasonable match to the system bandwidth and provides relatively good amplitude coverage. The next question we must ask is: how different is the output of the system when the input is sampled and held for five seconds from the output when the original signal is input. This corresponds to the problem of showing that the (A) to a loop in Figure 3.1 is the same as the (A) to (B) to (b) to (a) loop. One way to approach this problem is in the frequency domain. Figure 3.20 shows a sine wave of frequency 0.01 Hz., a representative input. The dashed line corresponds to the five second sample and hold signal. Figure 3.21 shows the difference between these two inputs and Figure 3.22 show the power spectral density of the difference signal. Note that all components except the fundamental at .01 Hz. would not be passed by the linear subsystem and so only the .01 Hz. signal contributes to any 121
PAGE 129
150 100 "' "' ;:;: 50 0 0 200 "' "' ;:;: :::> 100 z 0 0 150 100 150 200 AMPUTUDE BIN r ,h 20 40 60 AMPUTUDE BIN "' "' ;:;: :::> z "' "' ;:;: 200 100 0 0 200 :::> 100 z 0 0 ill l 20 40 60 AMPLITUDE BIN rnil 20 40 60 AMPLITUDE BIN Figure 3.19: Histogram of the input amplitudes for one second hold (upper left) and of the output amplitudes for one second hold (lower left), five second hold (upper right) and ten second hold (lower right). 122
PAGE 130
1 0.9 0.8 0.7 0.6 {/) E< 0.5 z ::> 0.4 0.3 0.2 0.1 0 0 ,, 200 400 600 TIME(SECS) 800 ,. ,, 1000 ,, 1200 Figure 3.20: Sine wave input of frequency 0.01 Hz. and the corresponding five second sample and hold. 123
PAGE 131
0.2 0.15 0.1 0.05 {/} ,.. z 0 :::> 0.05 0.1 0.15 0.2 0 200 400 600 BOO 1000 1200 TIME(SECS) Figure 3.21: Difference between original sine wave and its five second sample and hold. 124
PAGE 132
1 0.9 0.8 0.7 0.6 {/] ... 0.5 z :::> 0.4 0.3 0.2 0.1 0 I I 'II l 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 FREQ(HZ) Figure 3.22: PSD of difference signal. 125
PAGE 133
difference in the outputs. However, we have normalized the units in Figure 3.22 and if we carry through the calculations in absolute units we find that the ratio of the difference signal power to original sine wave power is .003 or 50 dB. We believe this is insignificant. We can also carry out the comparison in the time domain. Figure 3.23 shows the output of the simulated system when the sine wave sampled at 0.1 second is input (solid line) and when the five second hold signal is input. We do all processing digitally and so there is always a hold time involved but clearly 0.1 second is very similar to a continuous signal for this system. The two outputs are nearly identical again supporting the use of five second holding times. It is also interesting to note that as frequencies get higher the difference is greater. In fact, at very high relative frequencies we find the phenomenon shown in Figure 3.24. In that figure the top plate shows a sine wave sampled and held at a sample time too large, but in such a phase with the signal so that the samples occur at the extremes of the signal. The lower plate show the same signal and the same sampling time, but with the phase displaced. In this case the power of the sampled signal can be much lower than the true signal if we sample out of phase with the input. This is the time domain analogue of the picket fence effect.' 1 "Picket fence" is usually used to describe a frequency domain phenomenon. The discrete Fourier transform can be thought of as observing the Fourier transform with a set of bandpass filters. If the spectral component is not at the center frequency of a filter it is attenuated. This view through spaced bandpass filters is like observing a scene through a picket fence. 126
PAGE 134
"' ..... z ;:o ..... ;:o "" ..... ;:o 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 v 50 100 150 v v v 200 250 300 350 400 450 500 TIME(SECS) Figure 3.23: Output of the system when input is a sine wave (solid line) and when the input is the five second hold (dashed line). 127
PAGE 135
1 0.5 , ' ' , _j ' \ ' l ' j ' j "' ' t:: i i 0 \ \ z l i __ __j ::> i \ ___ _j 0.5 ____ j \ \ __ __j ____ j 1 0 100 200 300 400 500 600 700 800 900 1000 TIME(SECS) 1 0.5 "' \ ! \ ,_. ' _l 0 ' z \ ' \ _) \ \ J ::> _, ' ' 0.5 1 0 100 200 300 400 500 600 700 800 900 1000 TIME(SECS) Figure 3.24: Illustration of effect of sampling relatively high frequencies. 128
PAGE 136
However, since these relatively high frequency signals are greatly attenuated by the linear subsystem anyway, they contribute little to the output and can safely be ignored. The next question we must answer is how frequently to resample in going from (B) to (C) in Figure 3.1 and, since we sample both input and output syn chronously, how frequently to sample the output at (b). We claim that nothing is gained by resampling at a rate faster than one sample every five seconds because, as in the previous discussion, no new frequency information is gained when we view the overall system. It is true that by sampling the input at, say, one sample every second we can determine the presence of higher frequencies in the input than we would see by sampling at one sample every five seconds, but these higher frequencies are not passed by the linear system and so they are irrelevantthat is why we have chosen five second sampling times in the first place. We further note that nothing is gained by producing output samples more frequently than one every five seconds because, as discussed earlier in this chapter, the output is bandlimited and so can be reconstructed by interpolation from the five second samples. Also, such oversampling compromises the identification algorithm. One way to analyze this is to note that multiple sampling of the same input level produces multicollinearity. This is fully discussed in [6], but the effect is that sampling too frequently means that the samples are no longer independent and so the
PAGE 137
To summarize, we have developed a method to identify the system. We will use as input a constant switching pace input with sample and hold times of five seconds. We will resample the input at a rate of five seconds. For the identification algorithm we will sample the output at five seconds and construct a polynomial mapping between the discrete input and the discrete output. We will use two general formsone will consist of polynomials made from lagged inputs only (the Volterra series) and the other will consist of polynomials made from lagged inputs and outputs (the polynomial NARMAX model). We will choose the coefficients for the polynomial series so that they minimize the mean square output error. Once we have constructed the model it will operate by taking an arbitrary input and sampling it at five second intervals. These samples will be the input to the discrete model and will produce a discrete estimate of the output every five seconds. We will then construct the continuous output estimate by interpolation using the MAT LABTAI function file interp from the Signal Processing ToolboxTAI [28]. NARMAX Applied to the Example We will first consider the N ARM AX model applied to the above system. Before we apply our algorithm we present in the following paragraph an analytic derivation of a model for the system. This will both provide insight into the way systems work and provide another model which we can compare with the model the algorithm constructs. 130
PAGE 138
The linear subsystem we are dealing with is determined by the differential equation lOy+ y = x. If we discretize this into one second steps we obtain the discrete model y(t + 1)y(t) + 2y(t) = 2x(t). 1 10 10 Solving for y( t + 1) we get y(t + 1) = .9y(t) + .1x(t). Iterating this five times we get the analytic model for the linear portion to be y(t + 5) = .59y(t) + .41x(t). We then input one of the CSP test signals into the linear subsystem of the Wiener clipper system and produce an inputoutput data record. We then use this data record with the least squares algorithm (implemented using proceedure REG from SAS/ STAT6> 2 [38]) and with regressors y(t) and x(t) to get y(t + 5) = .585y(t) + .394x(t). This is the model of the linear subsystem constructed by our method and is in good agreement with the analytic result. This output from the linear subsystem then becomes the input to the nonlinear clipper. In order to determine the form of the output of the clipper we 2Registered trademark used to identify products of SAS Institute Inc. 131
PAGE 139
approximate its inputoutput function by a polynomial. Earlier in this chapter we noted that the nonlinear system produces responses at frequencies j, 2f and 4f when stimulated with a sine wave of frequency f and argued that this meant that its polynomial representation would be fourth order. (Recall Figure 3.11.) We also noted that the power at the successive coefficients decreased rapidly and argued that the coefficients of the successive powers in the power series representation must also decrease rapidly. The result of minimizing the square error between the clipper inputoutput function and the quadratic z ax+ bx2 over a set of points x uniformly distributed between zero and one is z = 1.31x.807x2 A plot of this quadratic with the clipper inputoutput function is shown in Fig ure 3.25. Since the nonlinearity is zero memory the output at time t resulting from input y( t) should be z(t) = 1.31y(t).807(y(t))2 If we substitute the form derived above for the output of the linear subsystem y(t + 5) = .59y(t) + .41x(t) into this we get that the overall system output at t + 5 should be z(t + 5) = 1.31(.59y(t) + .41x(t)).807(.59y(t) + .41x(tW Expanding this and collecting terms and then centering the variables about their 132
PAGE 140
0.6 0.5 0.4 "' 0.3 z ;::> ,.. ;::> 0. 0.2 ,.. ;::> 0 0.1 / 0' 0.1 0 0.2 0.3 0.4 0.5 0.6 0.7 0.9 1 INPUT UNITS Figure 3.25: Inputoutput function for the clipper and its quadratic approximation. 133
PAGE 141
means we get that In the true identification experiment we do not have access to the output of the linear subsystem and so do not know the y's. However, in the specialized environment we have set up we know them and it is interesting to regress the output z(t) on them and compare it to the above. The result is Clearly the two agree welL Note that we must use centered variables to discriminate between x and x2 and y and y2 Also note that the largest discrepancy is in the term because the true system has a contribution (although small) from y4 while the derived system includes only the quadratic approximation. Since y'1 is correlated with y2 and this is not removed by centering, the y4 effects are lumped in with the y2 terms. The above suggests that our idea of how the system works is correct; however, it doesn't help identify the system since we do not have access to the intermediate variables. We now we return to the problem of modeling the system as polynomials in factors we have access to, namely, the inputs and the outputs of the overall system. We first construct the data record by inputing a CSP signal with hold time of five seconds and amplitude values chosen from a uniform distribution. This produces an output. We then sample both the input and the output at the switching times (every five seconds) of the input signaL From this inputoutput 134
PAGE 142
record we construct the regressors. In this first example we will make them very simple and use only z(t), x(t) and x(t)z(t). Thus, the model equation is z(t + 5) = az(t) + bx(t) + cz(t)x(t). Our data record produces an observation of this every five seconds and we solve the resulting set of equations in the least squares sense as described before. We thus force the least square algorithm to build a one time step bilinear model [34]. We will discuss methods for selecting regressors from the larger set of Volterra or NARMAX polynomials more fully in Chapter Four. The result is the bilinear model z(t + 5) = .786z(t) + .487x(t).677z(t)x(t). All of the coefficients are included at a p level greater than .0001. At first this seems a very simple model, however, the R2 for it is .99069 and the Cp is 3.0. 3 All of this suggests that the model fits the original data well. Figure 3.26 shows the result of applying the model and the original system to another constant switching pace input with the same statistics as the training set. They clearly agree well and the standard deviation of the error is .057 units. In Figure 3.27 we compare the model and the original system for sine wave input. We see from the sine wave data that the agreement of the model and the underlying clipper system 3See Miller [32] for a full discussion of these. R2 is sometimes explained to be the portion of the variance in the data accounted for by the regression and Cr is an estimate of the cost of leaving regressors out. In "good" models it is near the number of regressors. Neither are completely reliable in this setting, but they both suggest that the model is valid. 135
PAGE 143
0.6 tJ] t: z ::> 0.2 .... N!fY/. v V .' v V ( .... In 'v \ 0 0 100 200 300 400 500 600 700 BOO 900 1000 TIME(SECS) Figure 3.26: Bilinear model (dashed line) and original system (solid line) applied to CSP input. is not perfect. However, the clipper is a much more "dramatic" nonlinearity than the MichalisMenton nonlinearities in our full simulation model. Also, the three term NARMAX model is very simple, yet it is clearly capturing properties of the underlying system. Volterra Series Applied to the Example Now we apply the method to the Volterra representation. We use the same inputoutput data as before only now we allow the regressors to be only lagged polynomials of the input with no lagged output terms. As above we use SASj ST AT
PAGE 144
UJ !'< z 0.5 ::0 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) UJ !'< 0.5 z ., ., , ., ' ' ::0 ' ' ' // /l ,/' ' ' ' ' ' ' I ' ._, 0 0 50 100 !50 200 250 300 350 400 450 500 T!ME(SECS) Figure 3.27: Model output (bottom platesolid) and system output (bottom platedashed) for sine wave input (upper plate). 137
PAGE 145
criteria. The discrete model produced is z(t) = .77x(t5) + .37x(t 10) + .17x(t 15).32x(t5)2 .07x(t10)2 .29x(t5)x(t10).16x(t5)x(t15). This is the Volterra polynomial model of the system. The R2 value is again above .99 and all coefficients have a p value above .001. There are two further observations we make about this result. Figure 3.28 shows the output of the model and the clipper system when the input is a CSP signal with the same statistics as the signal used to construct the model but with a different random number generator seed. The standard deviation of the error in this case is .069. There are two further observations we make about this result. One is that it involves seven terms instead of the three in the N ARMAX model yet the standard deviation of the error is higher in the Volterra model. This will be important in Chapter Four where the memory of the system is longer and the resulting number of coefficients grows accordingly. The other observation is that the quadratic terms included in the model weight crossterms more heavily than squares. For example, there is a crossterm between x(t5) and x(t10) with weight .29 and a crossterm between x(t5) and x(t15) with weight .16 while the square term x(t 10)2 has the low weight of .07 and the square term x( t 15 )2 does not appear at all. At first this seems odd, and to understand it we resort to an argument that uses intermediate variables as we did in the N ARM AX case. If we build a model of 138
PAGE 146
the linear subsystem using Volterra regressors we get y(t) = .43x(t5) + .24x(t10) + .14x(t15) + .07x(t20) + .04x(t25), where x is the input and y is the intermediate output of the linear subsystem. When we ptlt this in the equation for the nonlinearity N(z) = 1.31z.807z2 we see that the quadratic coefficients arise by multiplying corresponding coeffi cients in the linear model. Since the coefficients decrease rapidly and since,for example, .43 .24 > .242 and even .43 .14 > .242 we see why the crossterms predominate. Figure 3.29 shows a plot of the output of the model and the system under a sine input. The sine curve is interesting in that it shows the presence of the first harmonic as the ripple in the top part of the model output. This is expected since we allowed quadratic regressors and these produce the first harmonic as discussed earlier. Conclusion In conclusion, in this chapter we have constructed a method to build polyno mial models from inputoutput data and applied it to a nonlinear system. We have determined a class of inputs to use for determining the models and con structed a method for building the models from the resulting inputoutput data. 139
PAGE 147
0.6 [/) E< 0.4 z ::0 0.2 v 0 100 200 300 ' ' " 400 500 TIME(SECS) 600 700 800 900 1000 Figure 3.28: Volterra model output (dashed) and system output (solid) for CSP input. 1 '\ 1\ '\ f\ (\ 0.5 0 ,/ 200 400 f\ f\ \, v 600 TIME(SECS) \} \I 800 \ '\ \; \ \. 1000 1200 200 400 600 TIME(SECS) 800 1000 1200 Figure 3.29: Volterra model output (bottom platesolid) and system output (bottom platedashed) for a sine wave input (upper plate). 140
PAGE 148
The models that we have found have been quite simple but their agreement with the system output has been good, especially when we consider the constraints of limited time input and amplitude we have imposed. In Chapter Four we will apply these methods to the overall simulated system of Chapter One. 141
PAGE 149
CHAPTER 4 IDENTIFICATION OF THE MODEL Introduction In this chapter we will construct a model of the full simulated system presen ted in Chapter One. Clearly this simulated system is much more complicated than the systems considered in Chapter Three. In the full system there are three Wiener subsystems of the type considered previously that interact with each other. In Chapter Three we developed methods that worked relatively well in modeling these subsystems. A major addition to the full system is the venous pool. This venous subsystem receives inputs from the muscle subsystems and the visceral subsystems and acts as a flowdependent delay and as a mixing chamher. Functionally it acts as an integrator in that it averages the inputs with the previous state of the pool and smooths the output response. This subsystem introduces two new features that complicate the identification problem. One is that if we consider the input to the venous subsystem to be the oxygen concentration of the blood flowing into the chamber and the output to be the oxygen concentration of the blood flowing out of the chamber, we see that this system has hysteresis. This is illustrated in Figure 4.1. In that plot the lower leg of the loop represents slowly increasing the oxygen concentration of the input mixture 142
PAGE 150
into a chamber initially filled with blood of low oxygen concentration. Thus, the output concentration lags as the richer input slowly increases the concentration in the chamber. Then in the top part of the loop the reverse occurs: as the input mixture leans the output concentration decreases more slowly to yield the hys teresis loop. The output of this system depends not only on the current inputs to the system but also on the state of the mixing chamber. This state in turn depends on previous inputs and so, while it is possible to model such a system with lagged polynomials, it is clearly harder than the cases previously considered and may require more terms and more complicated terms. A second major difference between the full system and the systems of Chapter Three is the fact that the venous subsystem acts as a variable delay system where the transit time depends on other system variables. Again, this type of system can conceptually be modeled by lagged polynomials. For example, suppose that the system input is x(t), where 0 :S x(t) :S 1, and consider the system modeled by y(t) = x(t1)x(t2) + (1x(t1))x(t3). When x(t1) is near 1 (ie. high), the system output is nearly y(t) = x(t2) and when x(t1) is near 0 (ie. low) the system output is nearly y(t) = x(t3). 143
PAGE 151
120 100 80 60 40 20 0 12 3 4 56 7 8 9 10 INPUT Figure 4.1: Output from the mixing chamber vs input into the chamber showing hysteresis. 144
PAGE 152
This can be thought of as a variable delay system that clearly has a lagged polynomial representation. As in the case of the mixing chamber, this system can be represented as lagged polynomials, but again at the expense of using more complicated expressions. One important question we will consider in this chapter is how well the modeling method of Chapter Three extrapolates to these cases. The fact that we are using the simulation will let us examine the output of the subsystems directly and determine how well the modeling method is working and where it is failing. This would be almost impossible using real, noisy physiologic data. We begin this chapter by developing the NARMAX model. This involves evaluating the system in the same way as in Chapter Three by determining its input and output frequency properties and its amplitude properties. In Chapter Three we considered Wiener systems composed of cascades of linear subsystems and zero memory nonlinear subsystems. We saw that the test inputs must be designed to insure that the output of the linear subsystem has enough amplitude excursion to insure adequate testing of the following nonlinear subsystem. In the full system we have not only subsystems composed of Wiener systems but these are further cascaded with the venous pool and we must be careful to insure adequate amplitude excursions at each stage in the system. We then apply the modeling method developed in Chapter Three to determine N ARMAX models for the simulated system and discuss their validity. Next we develop a Volterra model for the simulated system. In addition to 145
PAGE 153
the frequency and amplitude information developed in the NARMAX section we will need to determine the approximate memory of the system and we will develop methods for doing this. We will then apply the methods and determine Volterra models for the system and discuss their validity and compare them to the NARMAX models. In the last section we will consider ways to construct models when the data contains noise. In Chapter Two we developed a method for modeling noise in the simulation and we will use this method to determine how the model building technique degrades in the presence of noise. We will then discuss methods for modifying the model building technique to include noisy data and present examples showing how the modification technique works. The NARMAX Model In this section we will build a NARMAX model of the simulated system. As m Chapter Three we will use the staircase like constant switching pace inputs with amplitude values drawn from a uniform distribution. In the case of the clipper nonlinearity in Chapter Three the amplitude limits were arbitrarily set to range from zero to one. Here the amplitude limits are determined by the physiologic parameters of the simulated subject. In a real laboratory experiment we ask the subject to warm up briefly and then increase the braking of the bike until the subject can barely pedal. This determines the upper limit for the load amplitude. In the classic V02max experiments, in which the load is constantly 146
PAGE 154
ramped upward, the subject quickly fatigues at this level and the experiment ends. However, under the dynamic inputs that we use, there is amplitude modulation imposed by the bicycle subsystem and so the subject does not see the entire amplitude excursion of the input. Also, the subject only experiences the higher loads for a brief time. Both of these effects combine to allow the subject to ride through this maximal amplitude. In the simulation, this is determined by the asymptotes of the inputoutput curves of the three Wiener subsystems discussed in Chapter Onethe input load to muscle oxygen consumption subsystem, the input load to muscle blood flow subsystem and the input load to cardiac output subsystem. In the examples in this chapter we will use the curves shown in Chapter One. Different experimental subjects would yield different curves and, of course, different values. In the present simulations the maximum load of the simulated subject is about 200 watts and we will use this as the upper input limit. The lower input limit is 0 watts. When we build the model we will center the input about the midrange input value (100 watts) and center the output about the output value when the input is at its midrange value. In the linear system case centering plays a less important role since there is no need to uncouple the effects of x and x2 as we discussed in Chapter Three. Secondly, in the linear case the mean of the output for an arbitrary signal is the output produced by the constant signal corresponding to the mean of the input signal. This is not true in the nonlinear case and explains why we center the output as described above. For this configuration of the 147
PAGE 155
simulation, the output is centered about 40.5 ml. oxygen per second. Billings and Voon [8] discuss this method 'of centering and presents examples. Now consider the hold times associated with the constant switching pace (CSP) inputs. This depends on the various time constants in the simulation. The realization we will study here is one in which all the time constants are 10 seconds. (In Chapter Two we discuss the effect of varying the time constants independently in the Step Inputs section.) According to Barstow and Mole [5] this is at the low range of possible values. Our data also suggest that most individuals will have larger time constants. However, the effect of larger time constants is just longer experimentsthey do not change the concepts involved. We feel that a parametric study of the effects of different time constants would be of value and we will suggest it as future work. Figures 4.2, 4.3 and 4.4 show plots of input in the time domain (upper left hand plate), output in the time domain (lower left hand plate), input in the frequency domain (upper right hand plate) and output in the frequency domain (lower right hand plate) for a one second hold time, a two second hold time and a five second hold time, respectively. The five second hold time meets the criteria discussed in Chapter Three of providing good bandwidth matching and good output amplitude excursion and we will choose this signal for the test input. As in Chapter Three we will choose the values at each time step from a uniform distribution with limits as discussed in the previous paragraphCSP input. Following the method developed in Chapter Three we input this CSP signal 148
PAGE 156
150 Ul !::: 100 ;z; ::> 50 1000 1500 TIME(SECS) 2000 60,.. 50 30 1000 1500 TIME(SECS) 2000 Ul !::: ;z; 1,, 0.5 FREQ(HZ) 1.. ::> 0.5 "' Ul 0.. orh 0 0.5 FREQ(HZ) Figure 4.2: Input level vs time( upper left plate), system output vs time(lower left plate), input PSD vs frequency (upper right plate) and output PSD vs frequency (lower right plate) for a hold time of one second. 149
PAGE 157
z 100 :::> 50 1000 1500 TIME(SECS) 2000 50 1000 1500 TIME(SECS) 2000 1,. 0.5 FREQ(HZ) 0.5 FREQ(HZ) Figure 4.3: Same as previous figure but for a 2 second hold time. 150
PAGE 158
150 100 ;:;> 50 u 1000 1500 TIME(SECS) 2000 1000 1500 TIME(SECS) 2000 0.5 FREQ(HZ) 0.5 0 0.5 FREQ(HZ) Figure 4.4: Same as previous figure but for a 5 second hold time. 151
PAGE 159
into the simulated system and sample both the input and the output at the switching times (every five seconds). A discrete model is then constructed by performing a least squares regression of the output on polynomials formed from lagged inputs and outputs. The next question to be answered is how to select the lagged polynomials to be used as regressors. We will first limit ourselves to linear and quadratic polynomials. As discussed in Chapter Three, this corresponds in the Wiener system case to approximating the inputoutput function of the zero memory nonlinearity with a Taylor series to squared terms. We will further select from the remaining polynomials by three methods. The first is called the FULL FITTED MODEL by SASjSTATffi [38] and in it the user specifies the regressors and the model is built using all of them. We will use this method in the following to build a very simple model consisting of bilinear terms only as we did in Chapter Three. The next method we will use is called STEPWISE by SASjSTATffi [38]. This method is sometimes called Efroymson's method and is discussed in Miller [32]. In short, this algorithm works iteratively by determining the remaining regressor that most reduces the residual sum of squares by its inclusion in the current model. It then calculates a ratio of this reduction in the residual sum of squares to the resulting residual sum of squares. If this value exceeds a threshold the regressor is included, if it does not it is not included. The method then determines the regressor in the model that least increases the residual sum of squares if it is dropped from the model. It then calculates a ratio of this increase 152
PAGE 160
in the residual sum of squares to the residual sum of squares with the regressor in the model and compares this ratio to another threshold. If the ratio is less than the threshold the regressor is dropped. This process continues until no new regressors are added and none dropped at any step. These ratios are called F values that we are dealing with the ratio of two chi square distributed random variables. We know that this is not the case, but as Miller [32] discusses, the method is of value and has theoretical justification even though the values are not normally distributed. It is useful in that it provides a stopping rule when selecting regressors. The final method we will use is called RSQUARE by SAS/ STAT [38] and consists of specifying a very large set of possible regressors and specifying how many terms we want in the final model. The method then considers all subsets of the specified size and chooses the one that minimizes the residual sum of squares. We will illustrate that this method can be used to find models that fit well but have many fewer terms than the models built by the stepwise method. There are many more possible methods of regressor subset selection but we will not use them here. The first model we will fit is a one time lag bilinear model: y(t) = ay(td)+ bx(t) + cy(td)x(t), where dis the hold time (5 seconds). The result of applying the least square technique to the data generated by the simulation as specified above and this 153
PAGE 161
model is that a = 42.7, b = 4.5 and c = 5.2, where we have normalized the input and output variables by theif approximate maximum excursions. The input variable varies from zero to 200 with center at 100 and so it is normalized by 100. The output variable varies from zero to about 80 and we will normalize it by dividing it by 50. The p value for each coefficient is .0001 and the R2 for the model is .9502. Thus, the resulting bilinear one step model is: y(t) = 42.7y(t5) + 4.5x(t)5.2y(t 5)x(t). ( 4.1) These values are then shifted by 40.5 since the output is centered as discussed above and then passed through the interpolation algorithm as discussed in Chap ter Three to produce output at the required time granularity. The resulting model is then tested by inputing various input sequences into it and into the simulation system from which it was built. Since we know that the model errors do not satisfy the assumptions used in deriving such statistics as Mallow's Cp and Akaike's information criterion, we do not use them at all in evaluating the model but instead rely on comparing model and underlying system outputs under the various test inputs. Again, recalling Figure 3.1, these experi ments are done by taking the test input and sampling it every five seconds and using this as the input to the discrete model. The resulting output corresponds to the output sampled every five seconds and this is then interpolated to produce the model output. The test signal is also input to the simulation (underlying system) and produces an output. These two outputs are then compared qual154
PAGE 162
itatively and quantitatively to assess the quality of the model. When the test input is a staircase signal with the same statistics as the signal used to construct this model but with a different seed for the random number generator (ie. a different realization of the process), the root mean square error (the square root of the mean square error) is 3.11 units. In this test we assume the sampling rate is one sample per second and we compare the two outputs each second. This result is not perfect, but we feel it is quite good in view of the very simple three term model. In Figure 4.5 we present the results of using a sine wave as a test signal and in Figure 4.6 we show the results of using a CSP test signal. In both cases the agreement is not perfect, but again the model is very simple. We see, as predicted in the introduction, that this simple model can capture the gross characteristics of the output, but fails to follow the subtle effects of the variable time delay and the mixing chamber. This is seen for the sine input because the model does not follow the elbow in the descending loop of the output. In the next case, instead of restricting the lagged polynomial regressors to the three used above we will allow virtually any regressor of the form x(td), x(t2d), y(td), y(t 2d), and quadratic products of these, where d is the sampling time (5 seconds). Thus, for example, x(t15) would be a candidate regressor as arey(t25), [x(t10)j2, [y(t5)]2and y(t30)x(t35). In order to get an idea of how far back in time we must go, Figure 4. 7 is a plot of the crosscovariance of the output with the lagged inputs. Since it is the relative values we are interested in, the plot is normalized to a maximum value 155
PAGE 163
100 TIME(SECS) / .. \, 40 / \ I 20 : /' //' ,/ / \ ' ' / ' \ ' ': / ,/ / ' ,/ / I 0 50 100 !50 200 250 300 350 400 450 500 TIME(SECS) Figure 4.5: Sine wave input and resulting output for the bilinear, one time step model. Solid=model/dashed=true system. 156
PAGE 164
200 150 (/) 100 z ::> (/) z ::> 50 0 100 50 40 ' 30 lr w II 150 200 ' \ ,, ( ,, \ u 250 300 ,., \ ,rl "\ ,_, l 350 400 TIME(SECS) '\ / '' '/ ', /'/\ \/'/\ 450 v 500 550 / ,, ' _/ I : \ \/ \j A( 600 100 150 200 250 300 350 400 450 500 550 600 TIME(SECS) Figure 4,6: CSP input (top) and resulting output for the bilinear, one time step modeL Solid=true system/ dashed=modeL 157
PAGE 165
0.8 ""' ::> 0. 0.6 :5 I ""' :::> 0.4 0. ""' ::> 0 "' 0 0.2 (.) >< 0 0 20 40 60 80 100 120 TIME LAG(SECS) Figure 4.7: Plot of the crosscorrelation of the system with lagged inputs vs lag. (Plot is normalized by the maximum yvalue.) of one. We will discuss this further in the section on Volterra models, but this crosscovariance relates closely to how regressors are entered into the model by the least squares regression algorithm. The value is nearly zero by sixty seconds suggesting that greater lag times are not needed. We feel that this is consistent with the physiology also in that system states more than sixty seconds old have little influence on the current state of the system we are considering. Cardiac output at rest is about five liters per minute and so the time to transit the three liter venous pool would be 36 seconds. Almost all of the tests are done at higher 158
PAGE 166
loads and so the corresponding transit time is much less. In fact, as we approach the 200 watt input level, the cardiac output approaches twenty liters per minute for a transit time of only nine seconds. Since the time constants for the systems that feed the venous pool in the current simulation are chosen to be 10 seconds, we feel that effects more than sixty seconds old are long forgotten. However, in constructing the discrete model we allowed lags up to 15d or 75 seconds. In this case we allow the regressors to be selected by the STEPWISE option in the S AS/ ST ATI'l> procedure REG [38] as discussed above. Again we use the CSP signals as test inputs to construct the model. This method with the threshold for inclusion and deletion both set to 0.15 (default) produced a model with eighteen terms shown in Table 4.1. The R2 for the model was 0.993 and the p values for all coefficients was .0001 except for the x(t10)x(t15) term (p=.0002) and the y(t 15) term (p=.004). The resulting model is complicated but agrees quite well with the simulated system for a variety of inputs. The root mean square error for another staircase signal with the same statistics as the input used to construct the model but a different seed is 1.28 units taken at a sampling time of one second. Figure 4.8 is a plot of the output from the underlying simulation system (top plate) and the model (bottom plate) for a staircase input with five second hold time and a different seed than the original signal used to construct the model. The agreement is very good. Figure 4.9 shows similar results for a one second hold CSP test input signal and Figure 4.10 shows the results for a ten second hold CSP input test signal. In both cases the agreement is again 159
PAGE 167
term coef. term coef. y( t5) 37.90 y( t5)x( t5) 1.70 y(t15) 1.54 y( t10)x( t10) 5.5 x(t) 4.66 y(t10)x(t15) 3.8 x( t5) 1.71 x(t)x(t) I 0.83 x(t10) 1.72 x( t )x( t5) 1.22 x(t15) .47 x( t )x( t10) 0.44 y( t5 )y( t5) 8.52 x( t5 )x( t10) 1.09 y( t5 )y( t10) 7.50 x( t5 )x( t15) 0.79 y(t5)x(t) 2.09 x( t10 )x( t15) 0.37 Table 4.1: Eighteen term NARMAX model. 160
PAGE 168
200 150 "' !:: 100 z :::> "' ... z :::> 50 0 100 60 50 40 30 20 100 lf w ll; u iU { 150 200 250 150 200 250 "300 350 400 TIME(SECS) 300 350 400 TIME(SECS) / 450 500 550 450 500 550 JL( / / 600 600 Figure 4.8: Output of simulation system (bottom platesolid) and the eighteen term model (bottom platedashed) for a five second hold staircase input (top plate). 161
PAGE 169
50 40 30 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 TIME(SECS) 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 TIME(SECS) Figure 4.9: Output of the simulation (top plate) and the eighteen term model (bottom plate) for a one second hold staircase input. very good. Figure 4.11 shows the results when the input is a sine wave. Here we see that the model (dashed line in the figure) is really trying to model the fine structure in the output and that the agreement is much better than for the three term bilinear model presenter earlier (Figure 4.5). Figure 4.12 shows the results when the input is a step function. Again we see much better matching than for the three term bilinear model (See Figure 4.6). The step presented is from 50 watts to 150 watts and although the model matches the gross output well, it does not seem to match the fine structure of the upslope as well as might be hoped. 162
PAGE 170
200 150 U) f< 100 z ;:> 50 0 I II 0 200 400 600 BOO 1000 1200 1400 1600 1BOO 2000 TIME(SECS) f I 20 i 0 200 400 600 BOO 1000 1200 1400 1600 1BOO 2000 TIME(SECS) Figure 4.10: Ten second hold CSP input (top plate) and corresponding model output (bottom platesolid line) and simulation output (bottom platedashed line). 163
PAGE 171
200 150 "' ... 100 z ::> 50 0 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) 60 //' /.... ' "' 40 I ... z ::> 20 0 0 50 100 150 200 250 300 350 400 450 500 TIME(SECS) Figure 4.11: Sine wave input (top plate) and simulated system output (bottom plate, solid line) and model output (bottom plate, dashed line). 164
PAGE 172
200 150 [/) c 100 z :::> 50 0 50 100 150 200 250 300 350 400 TIME(SECS) 60 [/) c 40 z :::> ff 20 / 0 50 100 150 200 250 300 350 400 TIME(SECS) Figure 4.12: Step input (top plate) and simulated system output (bottom plate, solid line) and eighteen term model output (bottom plate, dashed line). However, Figure 4.13 shows the results for a step function from 30 watts to 180 watts with the upslope section enlarged. Here it is clear that the model is trying to match the fine structure. By testing at several initial and final values of the step we find that at some values the model matches the fine structure well and at other values it averages as in Figure 4.12. This reflects the fact, as discussed in the introduction, that these second order effects and effects that depend on the remote past are difficult to model. There is just too much variation between these remote causes and the current effect for the least squares based algorithm 165
PAGE 173
"' 40 t: z p 20 / 0 50 100 150 200 TIME(SECS) 250 300 350 400 Figure 4.13: Enlarged response to a 30 watt to 180 watt step input for the simulation system (solid line) and the model (dashed line). to resolve them reliably and with small error. However, this model does show that it is possible to approximate even these second order effects. We now try to find a middle ground between the very simple threeterm model with its relatively poor performance and the complicated eighteen term model with its relatively good performance by nsing the selection method in SAS/ STAT [38] called RSQUARE that finds the minimum residual sum of squares over all regressor subsets of a given size. We first note that this process can be confusing. When there is a small number of regressors there is a great possibility that one is fitting the data and not the underlying system. This is manifested by large variations in the coefficients as we vary the number of regressors. For example, as we go from the best three term model to the best four term model, the coefficient of y(t5), which is included in both models moves from 42.5 to 56.2. Without testing the output of each model against the 166
PAGE 174
standard test inputs and comparing this to the output of the simulated system, it is hard to tell if the least squares algorithm is modeling the underlying system or the data. One guideline is model stability. For example, a model of the form y(t) = 1.1y(t5), where we are returning to normalized coefficients for ease of illustration (ie. 1.1 corresponds to 50* 1.1 = 55), can be rejected immediately because it is unstable. For pure autoregressive systems, stability is determined by the distribution of the roots of a polynomial associated with the z transform of the system. In the above case this polynomial is 11.1z which has the single root .909. This is inside the unit circle and so the model is unstable. In the case of the complicated model structures, stability becomes harder to evaluate because the z polynomial coefficients depend on the input x. For example, the model y(t) = .5y(t5) + .7y(t5)x(t) has z polynomial 1 .5z .7x(t)z, where we z transform only y and consider x(t) = c, a constant. (ie. we consider the various autoregressive models under constant input.) For the constant input x = .1 the z polynomial is 1 .57 z and this model is stable but for the constant input x = 1 the z polynomial is 11.2z and the model is unstable. This is further complicated by the presence of terms like y(t5)2 and y(t5)y(t10), but a quick stability check like the above can help in rejecting some of the models built by the best subset algorithm we are now considering. This method still involves trial and error. We first decide the maximum time lag to allow and from the inputoutput data construct the extended data set that includes a large number of possible regressors. We then choose the size of the model we want to build 167
PAGE 175
term coef. y( t5) 35.11 x( t) 4.69 x( t5) 1.20 x(t10) 1.90 y(t10)x(t10) 5.45 x( t )x( t5) 1.40 Table 4.2: Six term NARMAX modelsee text for details. (the number of regressors). The SAS/ proceedure REG [38] then goes through all subsets of the specified size and gives the ones that fit the data best. We then test these for stability and their performance against other data sets. In general, we find that as the number of regressors increases, the models fit better. However, by using this algorithm and the methods discussed above we have found a sixterm model defined in Table 4.2. All coefficients have p value .0001 and the R2 for the model is .990. Testing this model with a fivesecond staircase input with a seed different than the input yields a root mean square error of 1.59 units when measured at a one second sampling rate. This is plotted in Figure 4.14 and this sixterm model also gives quite good matches to sine and step input. Thus, we have developed NARMAX models for the simulated system using 168
PAGE 176
200 !50 z 100 ::> 50 0 100 lr tU u r !50 200 { Ali v 250 300 350 400 450 500 550 600 TIME(SECS) 50 ,, / '' _, \ ' J / /'' I\ I'> ''' J ,, 100 !50 200 250 300 350 400 450 500 550 600 TIME(SECS) Figure 4.14: Output of simulation (bottom platesolid) and the sixterm model (bottom platedashed) for CSP input (top plate). 169
PAGE 177
several regressor selection techniques. The last method has produced a model with only six terms but which matches the underlying system to within a root mean square error at each one second sample point of only 1.59 units. We feel that this is very good performance for such a simple modeL As noted earlier V02max as determined by a ramp input is frequently used in practical applications. When we input a ramp to the simulation the V02max is 52.18 ml/sec. When we input this ramp into the eighteen term model we get the estimate of V02max to be 50.69 ml/sec and when we input this ramp to the six term model we get 57.33 ml/sec. Both of these predictions are within ten percent of the simulation value and the eighteen term model is within 2.9 percent. Although care must be taken when extrapolating the model beyond its valid range, we have found human V02max values to vary about five percent from day to day and so this is a good result. The Velterra Model In this section we will build Volterra models for the simulated system using the method developed in Chapter Three. As discussed there, Volterra models allow only polynomials in the lagged input variables as regressors. We will use the same type of inputs that we used in building the NARMAX model in the previous section. As in Chapter Three we will sample the input and output at five second intervals and use this data to build a discrete model using least squares to regress the output on the Volterra polynomials of the input. We will then obtain the final model by interpolation. There are two fundamental questions that must 170
PAGE 178
1 0.9 0.8 0.7 "' 0.6 .... z ::> 0.5 Cl "' 0.. 0.4 0.3 0.2 0.1 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 FREQ(HZ) Figure 4.15: Power spectral density of the sine wave input. be answered in selecting the Volterra polynomials to use as regressors: what is the highest order polynomial needed and what is the largest lag time needed. We address the first question as we did in Chapter Three by inputing a sine wave into the simulated system. The power spectral density of the input (PSD) is shown in Figure 4.15. As we have discussed in Appendix A the output of the system under this input consists of the fundamental and the harmonics of the input and from the magnitude of the harmonics one can determine the relative contribution of the different powers in the Volterra model. Figure 4.16 shows the corresponding power spectral density of the system output under the sine wave input. There is a large contribution at the fundamental, as expected, and a small 171
PAGE 179
1 0.9 O.B 0.7 [/) 0.61t: z :::> 0.5 Q [/) 0.. 0.4 0.310.2 0.1 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 FREQ(HZ) Figure 4.16: PSD of the output of the simulation system under sine wave input. contribution at the first harmonic and a very small contribution from the seconi::l harmonic. Figure 4.17 shows the other terms in the PSD renormalized so that the first harmonic has value one. The purpose of this is to show that other terms are indeed present but that they are of very small magnitude. As discussed in Chapter Three the first harmonic corresponds to quadratic terms in the Volterra series the second harmonic corresponds to cubic terms in the Volterra series and so on for higher harmonics. Since the relative size of the contribution from higher harmonics is small we will build models using only terms up to quadratic. The next thing to determine is the largest lag time needed. We briefly dis172
PAGE 180
1 0.910.8 0.7 "' E:: 0.6fz ::> 0.5 "' Ul 0.. 0.4 0.3 0.2 0.1 0 )\ 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 FREQ(HZ) Figure 4.17: PSD of the output shown in the previous figure with the fundamental deleted and terms renormalized so that the first harmonic has value one. 173
PAGE 181
cussed this in the section on the N ARM AX model and obtained the result using crosscovariance that about sixty seconds was a sufficiently long time. Another way to view this is that this maximum lag time is really the memory of the systern and another way to determine it is to measure how long it takes the system to forget its past. In order to determine this we make two different runs of the simulation using inputs that are very different at first but then become the same. Thus, at the time when the inputs become the same the system has been prepared by the differing inputs to be in different states. We then observe how quickly the outputs from the two experiments become the same. Tills measures how quickly the system forgets its previous states and thus is an estimate of system memory. Figure 4.18 shows the sample inputs in the top plot. We first construct a random input. This is the solid line in the top plate of the figure. We then reflect this input about its value at 250 seconds and use this as the second input from the starting time until 250 seconds. This is the dashed line in the figure. After 250 seconds the inputs are the same. The bottom plate of the figure shows the outputs corresponding to the above inputs. They come together at about 310 seconds, again confirming that the system memory is about 60 seconds. From the above argument we feel that all posible regressors will be included if the maximum lag time is 75 seconds. As discussed before, we consider only linear and quadratic Volterra polynomials. Fitting this all inclusive model using the S AS/STAT REG proceedure [38] with STEPWISE selection produces the model in Table 4.3 with eleven linear 174
PAGE 182
150 200 250 300 350 400 TIME(SECS) 60 ... __ / ...... ,\ _. ...... ,_ .... _____ ... /'\ / __ /.',\\... I 40 ', ,...... __ ... "' ['< \ ..... "" z ;:> 20 0 100 150 200 250 300 350 400 TIME(SECS) Figure 4.18: Differing inputs (top plate) and their corresponding outputs (bottom plate). 175
PAGE 183
terms and 21 quadratic terms. This is a very complicated 32 term model. First note that the largest lag time the selection procedure chooses is fifty seconds again strengthening the hypothesis that the system memory is approximately sixty seconds. When we test this model against the system with the standard fivesecond CSP signal with a different seed, we find that the root mean square error is 2.01 units when calculated for a onesecond sampling time. This is worse than for the six term N ARMAX model. Also, Figure 4.19 shows a sine wave input in the top plate and the corresponding simulation system output (dashed) and model output (solid) in the bottom plate. The match is not bad but clearly not as good as for the N ARMAX model. There is little evidence that the model is picking up the fine structure. One reason for this is shown in Figure 4.20. If we look at the Volterra model under a sine wave input we see that the linear terms produce a single sine wave output at the same frequency as the input but with differing amplitude and phase (top plate of figure). Similarly, the quadratic terms produce a single sine wave output with the frequency of the first harmonic of the input (bottom plate of figure). These two terms are summed to yield the overall output and these two sine waves simply are not enough to build the fine structure in the output. It seems that this fine structure must be contained in the higher frequency terms in the PSD plot of Figure 4.12. In summary, we have constructed a Volterra model of the simulated system. Much as in Chapter Three, this model has required more terms than theNARMAX model and has 'not produced as good a model. 176
PAGE 184
term coef. term coef. x( t) 4.67 x(t5)x(t15) 0.35 x(t5) 4.76 x(t10)x(t10) 0.61 x(t10) 5.11 x(t10)x(t15) 1.09 x(t15) 3.97 x(t10)x(t20) 0.72 x(t 20) 2.78 x(tlO)x(t25) 0.43 x(t25) 1.17 x(tlO)x(t30) 0.36 x(t30) II 0.74 x(t15)x(t20) 0.61 x(t35) 0.44 x(t15)x(t25) 0.49 x(t40) 0.25 x(t15)x(t30) I 0.37 x(t45) 0.20 x(t15)x(t35) 0.27 x(t)x(t) 0 .68 x(t20)x(t25) 0.37 x(t)x(t5) 1 1.37 x(t20)x(t30) 1 0.16 x(t)x(t10) 0.55 x(t20)x(t35) 0.41 x(t5)x(t5) 0.40 x(t25)x(t25) 0.39 x(t5)x(t10) 1.38 x(t25)x(t 30) 10.31 Table 4.3: Thirty twoterm Volterra model. 177
PAGE 185
UJ 150 3 100 50 200 300 400 500 600 700 800 900 1000 TIME(SECS) 100 200 300 400 500 600 700 800 900 1000 TIME(SECS) Figure 4.19: Sine wave input to the system and model (top plate) and resulting output for the system (dashed) and model (solid) (bottom plate). 178
PAGE 186
u 60 "' "' "'"" ::>1 40 f:;' :::> 0.. E< :::> 20 0 50 100 150 200 TIME(SECS) 0 u n "' 2 "' "'""' ::>1 4 f:;' :::> 0.. 6 E< :::> 0 8 v v 50 100 150 200 TIME(SECS) Figure 4.20: Decomposition of the Volterra model into the fundamental and first harmonic. 179
PAGE 187
Building Models in the Presence of Noise In determing the coefficients in the model so far we have considered only the noisefree case. In the last section of this thesis, we will consider extending the model to the noisy simulation as developed in Chapter Two. This is a difficult subject and the theory is not complete, but the simulation approach of this thesis is very useful in that we can see the results of modeling techniques for both noise free and noisy data. We will first determine the effect of introducting noise on the models. We will next discuss a method for correcting the models built with noisy data and study its efficacy with some specific examples. Finally we will discuss other methods for building models in the presence of noise. First consider the bilinear onestep N ARM AX model built in a previous sec tion of this chapter. Using noisefree data this model had coefficients as shown in equation 4.1. For this model the root mean square error when the model was tested against a five second hold CSP signal was 3.11 units. We now add Gaussian white noise of zero mean and variance 2.6 to the simulation output as discussed in Chapter Two. The bilinear onestep model constructed using this noisy data set has the coefficients shown in Table 4.4. The root mean square error in this case has gone up to 3.40 units. As expected, the modeling method degrades in the presence of noise. Now consider the sixterm NARMAX model built with the SAS/STATI& [38] proceedure REG and the RSQUARE method of model selection. The model 180
PAGE 188
term coef. y(t5) 33.80 x(t) 4.96 y(t5)x(t5) 1.40 Table 4.4: The three term bilinear NARMAX modelbuilt using noisy data. constructed using noise free data had coefficients as shown in Table 4.1. The corresponding model using noisy data as described above has coefficients as shown in Table 4.5. In this case the noisefree root mean square error was 1.28 units and the root mean square error with noisy data is 3.13 units. Again the performance of the modeling method degrades in the noisy data case. Now we consider a method for modifying the model constructed with noisy data. In the noisefree case we constructed the model: y(t) = qnl + e(t), (4.2) where y(t) is the observed output at timet, 0 is the column parameter vector,
PAGE 189
term coef. y(t5) 38.27 x( t) 5.13 x(t5) 0.73 x(t10) 1.73 y(t10)x(t10) 5.76 x(t)x(t5) 1.79 Table 4.5: The sixterm NARMAX modelbuilt using noisy data. to inaccuracies in the complicated measuring system that is used to determine breath to breath oxygen consumption measurement noise. If we view the system in a state space context we can write: X(t + 1) = F{X(t),x(t),e1(t)} y(t) G{X(t), e2(t)}, ( 4.3) ( 4.4) where X is the underlying state vector, x is the input to the system, F is the state evolution function, G is the observation equation, e1 is a vector of process noise and e2 is the measurement noise. In reality, both e1 and e2 are very hard to determine. For example, one component of e1 could be the noise associated with muscle oxygen consumption. As discussed in Chapter. One experiments to 182
PAGE 190
determine this are very limited. However, as discussed in Chapter Two, the subsequent cascades of the entire system are smoothing and such noise probably has little effect on the actual V02 output. In fact, the only real noise data we can measure without great difficulty is the overall output noise and, as discussed in Chapter Two, this can be approximated as additive, white, Gaussian noise. This is the approach we now consider and we assume that we can, at least conceptually, solve equation 4.3 for X(t) in terms of previous x and e1 values as X(t) = H{x,e1}. Putting this into the second equation ( 4.4) gives y(t) = G{H{x,e1},e2(t)}. We now make the assumption that the effects of e1 and e2 can be lumped into one additive noise term, e(t). This assumption is commonly made [8] and is reasonable here because the effects of noise early in the system are smoothed and the late effects, such as noise in cardiac output, are, to a first approximation, white and Gaussian as discussed in Chapter Two. This, when added to the presumed white, Gaussian effect of V02 measurement, gives an overall additive, white Gaussian noise term and allows us to write the inputoutput equation as y(t) = K{x} + e(t), where K is the appropriate composition of G and If above. This suggests that the model of equation 4.1 is still valid, only now e(t) includes process n01se, 183
PAGE 191
measurement noise and misfit. Now consider a model of the form in equation 4.2. As shown in [41] the least squares estimate of e, {J takes the form {J = t
PAGE 192
a ARMAX model as shown in the following example. Suppose the underlying system is given by y(t) = ay(t1) + by(t2). We actually measure a process z(t), which is an observation of y(t) contaminated by zero mean white Gaussian noise. Thus, z(t) = y(t) + e(t), or, y(t) = z(t)e(t). Thus, the observations satisfy z(t)e(t) = a(z(t1)e(t1)) + b(z(t2)e(t2)). Solving the above for z we get z(t) = az(t1) + bz(t2)ae(t1)be(t2) + e(t). If we lump the noise containing terms together as e(t), we have the least squares equation z(t) = az(t1) + bz(t2) + e(t), with two parameters a and b. However, the regressors z(t1) and z(t2) are clearly correlated with e( t) and the estimates will be biased. However, if we include e(t1) and e(t2) as regressors then the residual is e(t). This is a white Gaussian noise and so is uncorrelated with the regressors. This is the basis 185
PAGE 193
for the extended least squares method as presented in [2]. A similar argument holds for nonlinear systems [9], orily here there can be crossproducts of noise and regressors. For, example, suppose one of the regressors is [y( t 1 )J2x( t 1 ). Then substitution of y(t 1) = z(t 1)e(t 1) would give terms like z(t 1)e(t 1) and e( t 1 )2 However, we can include these terms as explicit regressors just as we included e(t1) and e(t2) in the linear case. We now apply extended least squares to the simulation system with nmse. The most common way to implement extended least squares is in the context of recursive least squares. In this method [15], an explicit process is established to update the least squares solution as each new data point is observed. This can be done either by a manipulation of the matrix that is used to calculate the least squares estimate or by applying the Kalman filtering algorithm. In either case, this allows one to calculate the parameter vector at time t, then use the inputoutput data available until time t to calculate the estimated output at time t + 1. This is then compared to the observed output at time t + 1 to obtain the current estimate of e. This method is applied iteratively to construct estimates of all e until the current time and they are then used as regressors as discussed above. We prefer to implement the algorithm in batch mode. With this approach, we first generate an estimate of (} by running the least squares algoritm on the noisy data using only the lagged N ARM AX regressors and no noise regressors. We then calculate the estimated outputs y(t) using this iJ and the corresponding 186
PAGE 194
residuals (y(t)y(t)). This becomes the estimate of e(t) and we use these to calculate the noise containing regressors for the first iteration. Least squares is then repeated using the full set of noise free and noise containing regressors and this gives us a first corrected estimate of B. This corrected value of {; is then used to calculate a second corrected estimate of y(t). A second corrected estimate of e can then be made and the process continues iteratively. We now apply this algorithm to the two cases presented earlier in this section. They illustrate the variety of results. In the first case, consider the three term bilinear model. In this case we know that there is a relatively large misfit error because the model has so few terms. When we apply the correction algorithm it gives after one step the model y(t) = 1.14y(t5) + .05x(t).002y(t5)x(5). This model is unstable as can be seen by the methods discussed earlier in this chapter. The problem is that the model diverges when the test input is applied as shown in Figure 4.21. This means that we can not calculate the residuals and the iteration stopsthere is no improvement in the model. On the other hand, if we apply the extended least squares method to the sixterm model we find that the method converges after nineteen steps to the corrected model shown in Table 4.6. Comparing this model with the noisefree model shown in Table 4.4 we see that they are quite similar. In fact, the model of Table 4 .. 6 fits the test signal used to produce it better than the noisefree model. We have found that 187
PAGE 195
0 100 ;:, 200 0 50 100 150 200 250 TIME{SECS) Figure 4.21: The output of the first corrected model in the three term bilinear case. term I coef. y(t5) 35.00 x( t) 4.99 x(t5) 1.10 x(t10) 1.72 y(tlO)x(t10) 5.90 x(t)x(t5) 1.70 Table 4.6: Sixterm corrected model. 188
PAGE 196
"' ""' 0.5 iii Figure 4.22: Distance (ZZ) between noise free e and successive extended least squares e (normalized) vs. iteration number. this can indicate another problem because the nmse correction algorithm can converge to a model that fits the test data well but does not fit other inputs. In this example, this is not the case and the corrected model gives good agreement with the underlying simulation under a wide variety of inputs. It is interesting to. see how the extended least estimates converge. The [2 distance between the noise free estimate of fJ and the estimate at each successive iteration relates to the bias of the noisy model. Figure 4.22 shows a plot of this 12 distance between the noise free estimate of table 4.4 and the evolving model vs. iteration number. (Recall that the values in the tables are multiplied by a scale factor of 50. In Figure 4.22 we have normalized the distance by its maximum value.) As the above results indicate the extended least squares algorithm can be very helpful in correcting models built from noisy data. There are many other methods that can be used. For example, modified maximum likelihood methods and modified 189
PAGE 197
instrumental variable methods have been discussed in the literature. The author is also developing a method that "first smooths the raw data first and then uses either regular least squares methods or any of the above methods. However, a full development and comparison of these methods will be left for future work. 190
PAGE 198
CHAPTER 5 SUMMARY AND FUTURE WORK Summary of the Thesis In this thesis we have considered the problem of building mathematical models of systems from inputoutput data. The motivation for this project has been the goal of finding a method for identifying the system of an exercising human where input is work load and output is the rate of oxygen consumption. This real system is difficult and expensive to study and so a simulation model based on physiologic principles and existing data was constructed. This model was presented in Chapter One. It was important to include all subsystems in this model and so models were developed for all parts of the experimental systemthe computer, the bicycle, the measuring devices and the subject. The subject subsystem was further modeled by defining several subsystems from basic physiologythe muscle system, the cardiac system, the lung system, the visceral system, the venous system and the arterial system. The muscle system and the cardiac systems were further modeled in terms of the Wiener model for nonlinear systems. The venous system was further modeled as a flowdependent buffer system. In Chapter Two the behavior of the simulation was determined under severa! inputs often used in systems analysis, such as steps, ramps, sme waves 191
PAGE 199
and random inputs. We compared the output of the simulation qualitatively to known outputs such as steps and quantitatively to experimental data. In all cases good agreement was obtained. In Chapter Two a noise model was also added to the simulation. Based on experimental observations we concluded that additive, white Gaussian noise was a realistic approximate noise model. We then turned to the problem of constructing an inputoutput model of a system such as the above. It was decided to consider only lagged polynomial models such as Volterra series and NARMAX models. In Chapter Three an identification method was constructed by sampling the input and the corresponding output and then using least squares regression to determine the coefficients of the lagged polynomials. We also considered a class of constant switching pace signals and discussed their benefits as input test signals. It was also necessary to develop a strategy for choosing sampling rates based on the bandwidth of the input and the output of the system to be identified. Finally in Chapter Three Volterra models and NARMAX models were devised for a simple Wiener system consisting of a first order linear system in series with a clipper nonlinearity. It was concluded that the N ARM AX model generally agreed better with the simulation result than the Volterra model and had the added advantage of having fewer terms. Chapter Four was devoted to applying the methods of Chapter Three to the full simulated system. We applied the least square algorithm using several criteria for deciding if a given lagged polynomial term should be included in the model. 192
PAGE 200
With these techniques several models of the simulation were constructed which showed good agreement between the model and the simulation. In Chapter Four we also developed a method, extended least squares, for constructing models in the presence of noise. It was shown that in some cases it yields good results and converges to good models while in other cases it does not converge at all. Overall, we have developed a method that will produce a lagged polynomial model of a system from inputoutput data. The method is applicable to the identification problem of the human work loadoxygen consumption system. In particular, the N ARM AX model selection that builds the best model of a given size can produce accurate and parsemonious models. Future Work The problem of developing a reliable, practical algorithm to identify the work loadoxygen consumption system of humans has proved to be large and difficult. It requires knowledge from diverse fields such as mathematics, physiology and electrical engineering. While writing this thesis the author has learned a lot about each of these subjects and has developed a feasible method applicable to the original problem; but there is much more to be done. We now list some of this future work. The first area of future work involves data collection. During the early part of this thesis we had access to a data. collection system designed for research. The input and output was easy to manipulate and sample and the bicycle itself 193
PAGE 201
was easy to control in a predictable manner. When one of the author's thesis advisors left the University of Colorado, this laboratory was closed and we were forced to continue research on commercial systems designed to determine V02max Although these systems work well at this specific task, they simply do not have the versatility re.quired for further work without modification. These systems control the bicycle through an A/D card that is controllable from the C programming language. The author is currently developing a data processing package that will allow much better control over input and data collection and this package will be used this to further validate the model, to determine parameters of the simulation and to collect data to test the method of system identification developed in this thesis. Regarding the simulation itself, as we have pointed out, data about each of the subsystems is very hard to obtain. It requires intervention into relatively inaccessible parts of human subjects during exercise. In spite of these difficulties there is a trickle of data about these subsystems and we will continue to use it to evaluate and modify the simulation and the simulation noise model as it becomes available. The identification method is also a source for future work. As discussed in Chapter Three some pseudorandom sequences have better properties than others for use in identification. The same will be true of the constant switching pace signals and we plan to develop a method to determine these best sequences. Another area of future research is the determination of a method of selecting regressors 194
PAGE 202
m the least square algorithm that works best. In this thesis we have applied these methods manually, but the process could be automated to find true best sets of regressors. We would also like to do study how the identification method performs under different simulation parameters as discussed in Chapter Four. Also, the area of developing models in the presence of noisy data needs further research. There are many methods suggested in the literature and as mentioned in Chapter Four, the results can be quite varied and there is no unified theory to predict applicability of the methods. We feel that the simulation provides an excellent tool to evaluate these methods and develop new ones. Finally, the relationship between the parameters in the simulation or the underlying human system and the coefficients of the lagged polynomials in the derived models is another area for future work. Just as switching from the time domain to the frequency domain by Fourier transformation gives a new perspec tive from which to view a signal, so switching from the inputoutput data to the lagged polynomial representation gives a different view of the system. The simulation is a very good tool for determining relations between these two representations and we hope to investigate it further. 195
PAGE 203
APPENDIX A Wiener vs Hammerstein Models In this appendix we will analyze the observations in the section titled The Muscle Subsystem about the input of sinusoids into the Wiener and Hammerstein models of the system. As shown in Figures 1.9 and 1.10, a single sine input with high frequency relative to the time constant of the linear subsystem yields a sinelike output in both models, but the Hammerstein model has lower mean output. On the other hand the input of a sine of low frequency relative to the time constant yields a periodic, nonsinusoidal output with mean again lower in the Hammerstein case. It is known [7] that a causal, linear system can be represented as a convolution: y(t) = f"' h(r)x(tr)dr where x is the input, y is the output and h is a function called the impulse response of the system. If we let x(t) = Aexp(jwt), we find that: y( t) = ('"' h( r )A exp(jw( t r) )dr = exp(jwt )A r= h( r) exp(jwr )dr Jo Jn The integral term is a constant and so a sinusoidal input to a linear system yields 196
PAGE 204
a sinusoid of the same frequency as output. (Note that the integral may be complex and, thus, the output can also involve a phase shift.) In the case of a nonlinear system there is a corresponding multiterm representation called the WienerVolterra series. (See Rugh [36] or Schetzen [39] for a full When this series has a finite number of terms, a converse of the above result holds: If the output resulting from all possible single sinusoids is a sinusoid of the same frequency as the input, then the system is linear. To see this write the WienerVolterra representation of the system y(t) =I h1h)x(tr1) +I I h2(r 1,r2)x(trt)x(tr2)dr1dr2 + +I I J hnh, Tz, .. Tn)x(tT1)x(t Tz) .. x(tTn)dr,dr2 drn where again x is the input, y is the output the h 1 are called the Wiener kernels, and the lower limit of integration is zero and the upper limit is Mthe system memory. Let x( t) = A exp(jwt) be a sinusoid of frequency w and under the hypothesis that the output is a sinusoid of frequency w and amplitude A we get Aexp(jwt) =I h1(r1)Aexp(jw(t ri))dr1 + I I hz(r,,rz)A2exp(jw((tr1 ) + (tr2))dr1dr2 + + l .. lhn(r,, .. ,rn)Anexp(jw((trl)+ .. +(t Tn))drl'"Tn = (exp(jwt))(A I h1(r1)exp(jwri)dr1 ) + (exp(2jwt))(A2 I hz(r1,r2)exp( jw(r1 + Tz)) + 197
PAGE 205
+( exp(njwt))(An J J hn( 1'1," Tn) exp(jw(TJ + .. + Tn) )dr1 .. drn) Since { exp( njwt)} is an orthogonal sequence for each interval [2k7r, (2( k + 1 )7r) we can show by an orthogonality argument that each integral term except the first is zero and since this holds for all frequencies, each of the higher order kernels is identically zero and so the system is linear. This knowledge and the results of testing only with high frequency sinusoids would lead one to believe that the systems we have developed are linear. Testing with sinusoids is very common in system analysis and the above analysis shows how important the choice of input can be in the design of a testing proceedure. Now consider the test with sinusoids in the frequency domain. Recall that a linear, first order subsystem is governed by ri; + y = i. In the Laplace domain this has the form r(sY(s)y(O)) + Y(s) = I(s) where capital letters represent Laplace transforms. Assuming y(O) = 0 we get H(s) = Y(s) = 1 I(s) rs+l where H( s) is called the transfer function of the system. By limiting ourselves to the frequency axis s = jw and multiplying H by its complex conjugate we get the power transfer function, which relates the power output at a given frequency 198
PAGE 206
1 0.9 0.8 0.7 "' 0.6 t: z :;, "' 0.5 "' "' 0 0.4 '" 0.3 0.2 0.1 0 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY Figure 0.1: Plot of power attenuation verses frequency w(radians/sec.) forT= 10. to the power input at that frequency: 1 P(w) = H(jw)"H(jw) = T 2 W 2 + 1 Figure .1 shows a plot of this function for T = 10 and Figure .2 shows the plot when T = 30. The effect of the time constant is that as it becomes larger the plot becomes narrower. This is often interpreted in terms of the bandwidth of the system which is defined as the frequency at which the power is one half the power at the center frequencyin this case w = 0. Figure .3 shows the power transfer function with the power spectral density (PSD )1 of a sinusoid (two spikes) 1 PSD is a measure of how the power of a signal is distributed in the frequency domain. See Benda! and Piersol [7] for a discussion. 199
PAGE 207
1 II 0.9 0.8 0.7 Ul 0.6 t: z :0 "' 0.5 "' "' 0 0.4 0. 0.3 0.2 0.1 1 0 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY Figure 0.2: Plot of power attenuation verses frequency w(rad.ians/sec.) forT = 30. 200
PAGE 208
LINEAR SYSTEM TRANSFER FUNCTION With "relatively high frequency sine 1.2,, p 0 w 0.8 e r u n i t s 0.6 1 0.5 0 Frequency(Hz) +Attenuation 0.5 Figure 0.3: Power transfer function discussed in text with the power spectral density of a sinusoid superimposed. superimposed and defines what we mean by "relatively high" or "relatively low" with respect to the time constant of the time constant of the systemr. The PSD of the output of a system can be determined by multiplying the PSD of the input by the power transfer function. In systems with high time constants (slowly responding systems) the bandwidth is low and the "relatively high" frequencies are attenuated while the "relatively low" frequencies are not. With this perspective consider the effect of the simulation on input sinusoids. Figure .4 shows the frequency domain representation of a single sinusoid of frequency w that is "relatively high" with respect to the time constant of the linear 201
PAGE 209
1 0.9 0.8 / \ / \ / \ ' ' I \ ,' \ 0.7 "' 0.6 !:: / \ z ;:> "' 0.5 "' "' 0 0.4 "" \ ' ' / \, I \ I \ 0.3 //' \\, 0.2 0.1 /' ',, / __ J_ __ 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY Figure 0.4: PSD of a sinusoid with "relatively high" frequency superimposed on the power transfer function discussed in the text. system. The power spectral density of the sinusoid is represe11ted by the spikes and the power transfer function of the linear system is superimposed. Now consider in the frequency domain the results of inputing this "relatively high" frequency sinusoid into a Wiener system and then into a Hammerstein system. In the Wiener case, shown again in Figure .5 the sinusoid will first enter the linear system and the output will be a sinusoid at the same frequency as the input, attenuated by about 0.2, as can be seen from Figure .4. The input in the time domain is shown in Figure .6. Let i(t) = Aexp(jwt) + M, 202
PAGE 210
LS io(t) I NLS L______ o(t) Figure 0.5: Wiener system block diagram. INPUT VS TIME 180 160 140 "' .... z ::> 120 .... 0. 1:: 100 80 60 0 50 100 !50 200 250 300 TIME(SECS) Figure 0.6: Sinusoidal input into the Wiener system in the time domain. 203
PAGE 211
OUTPUT OF THE LINEAR SYSTEM VS TIME 135 130 125 120 I 0 50 100 150 200 250 300 TIME(SECS) Figure 0. 7: Output of linear system part of the Wiener model. whereA is the amplitude and M is the mean value of the sinusiod. Then, using the notation of Figure .5, the output of the linear system is io(t) = Aexp(jwt + cp) + M, where A is the amplitude, cp is the phase, and M is the mean value of the output. This output is shown in Figure 7. Now let the zero memory nonlinear subsystem of the Wiener system be represented by N(). This means that if io(t) is input into it, then the output is N(io(t)) and assume that N(t) can be fit with a polynomial series of degree n. This assumption holds for the continuous curve nonlinearities we discussed in The 204
PAGE 212
Subject section of Chapter One and can be applied to discontinuous nonlinearities such as thresholds and dead zone by smoothing or by developing the series in terms of orthogonal polynomials such as Laguerre or Hermite polynomials where the coefficients are determined by integration. We will instead use a Taylor series applied to the MichaelisMenton approximation discussed in The Subject section: Since this function is smooth in the region we are interested in this local technique is valid. In fact, the only singularity occurs at k2 We need our expansion to be valid from near zero to near 200 watts and since k2 is positive we can expand about 100 and it will be valid in the desired region. If we do this expansion we get, after some manipulation z 100 z 100 2 k2+100(k2+100) + .. ). Now assume that M = 100. (The following works for other values of M, but the algebra is much more difficult.) The input into the nonlinear system (output from the linear system is z(t) = io(t) = 100 + Aexp(jwt + cp). The corresponding output from the nonlinear system, which is the overall output of the Wiener system is o( t) = k k11 00 3 + k, + 100 k1 k2 )(A exp(jwt + 'P (A exp(jwt + 'P) )' + .... k, + 100 k, + 100 k, + 100 205
PAGE 213
1 0.9 O.B 0.7 U1 0.6 t: z "' "' 0.5 "' "' 0 0.4 0. 0.3 0.2 0.1 I 0 I 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 O.B 1 FREQUENCY Figure 0.8: Frequency domain output of Wiener system"high frequency" case. Fundamental magnatude normalized to one. 207
PAGE 214
OUTPUT VS TIME FOR WIENER SYSTEM 44 43.5 43 "' 42.5 Ez :::> E42 :::> 0. E:::> 41.5 0 41 40.5 40 0 50 100 150 200 250 300 TIME(SECS) Figure 0.9: Time domain output of Wiener system"high frequency" case. 1 io(t) I Ht> __ L_S ___j o(t) Figure 0.10: Hammerstein model. 208
PAGE 215
Ul E< z ::> "' "' i>' 0 0. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 I I 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 FREQUENCY Figure 0.11: Output of the nonlinear part of the Hammerstein model. Fundamental amplitude normalized to one. fundamental at the frequency of the input sinusoid and the corresponding harmonies. Again, the amplitude drops off quickly as discussed above but harmonic components still remain. This can be seen in Figure .12, which shows the output of the nonlinear part of the Hammerstein system in the time domain. It is clearly not sinusoidal and still has harmonics present. This signal then becomes the input to the linear system in the Hammerstein model. As before we determine the output of the linear system by superimposing its power transfer function on this input. We do this in Figure .13 and see that the harmonics are attenuated even further. The time domain output is shown in Figure .14. 209
PAGE 216
OUTPUT OF NONLINEAR SYSTEM VS TIME 46 44 42 40 "' 38 Ez ;::> 36 34 32 30 28 0 50 100 150 200 250 300 TIME(SECS) Figure 0.12: Effect of the nonlinear system part of the Hammerstein model. 210
PAGE 217
UNEAR SYSTEM: TAU=10 1 0.9 0.6 0.7 UJ 0.6 ... z :::> "' .0.5 "' "' 0 0.4 0.. 0.3 0.2 0.1 0 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY(Hz) Figure 0.13: Power transfer function for the linear system superimposed on the power spectral density of the output from the nonlinearity in the Hammerstein system. 211
PAGE 218
OUTPUT VS TIME FOR HAMMERSTEIN SYSTEM 44 43 42 Ul z :::> ,.. 41 :::> 0. ,.. :::> 0 40 I I I 39 38 0 50 100 150 200 250 300 T!ME(SECS) Figure 0.14: Time domain output of Hammerstein model"high frequency" case. 212
PAGE 219
1 0.9 /\\ 0.8 / \ ' ' ' I I l \ ' ,' \ 0.7 ' / \ Ul 0.6 t: z :::> "' 0.5 t"' loo 0 0.4 p.. 0.3 0.2 0.1 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY Figure 0.15: Frequency domain representation of "low frequency" sinusoid. Thus, in both the Wiener and the Hammerstein cases for these "high frequency" inputs we have little of the harmonics remaining and the output is nearly sinusoidal in both cases. Also, in the Wiener case the harmonics are attenuated only by the decrease in the Taylor coefficients, while in the Hammerstein case they are attenuated by both the coefficient decrease and the effect of the linear filter. This explains why the relative amplitude is lower in the Hammerstein model case. Now consider a "relatively low" frequency sinusoid as shown in the frequency domain in Figure .15 and in the time domain in Figure .16. In the Wiener 213
PAGE 220
INPUT VS TIME 180 160 140 UJ ,... z ::> 120 ,... \ \ ::> I 0. ::: I 100 v 80 60 0 50 100 150 200 250 300 TIME(SECS) Figure 0.16: Time domain representation of "low frequency" sinusoid. 214
PAGE 221
OUTPUT OF THE LINEAR SYSTEM VS TIME 160 !50 140 130 "' t: z :;, 120 110 I v \ \ \ 90 0 50 100 !50 200 250 300 TIME(SECS) Figure 0.17: Output of the linear system for the Weiner modellow frequency case. model the linear system acts first and produces a relatively small attenuation as can be seen in Figure .15. The time domain output of the linear subsystem is shown in Figure .17. The nonlinear system then acts in a manner consistent with the coefficients of the Taylor expansion. Since the linear system attenuates the fundamental only slightly in this "low frequency" case some harmonics still have enough power to contribute noticibly to the output and we get a clearly nonlinear output as shown in Figure .18. In the Hammerstein case the nonlinear system acts first producing the time do215
PAGE 222
OUTPUT VS TIME FOR WIENER SYSTEM 42 "' t:: z ::> ,.. 40 :::> 0.. ,.. ::> 0 38 36 0 50 100 150 200 250 300 TIME(SECS) Figure 0.18: Output of the Wiener system"low frequency" case. 216
PAGE 223
OUTPUT OF NONLINEAR SYSTEM VS TIME 46 I\ I\ ( \ I \ (\ I \ 44 42 40 38 36 34 \ 32 30 v v v v v 0 50 100 150 200 250 300 TIME(SECS) Figure 0.19: Time domain output of the Hammerstein nonlinear subsystem"low frequency" case. main output shown in Figure .19 and the frequency domain output of Figure .20, where we have again normalized the fundamental magnatude to one. However, when the linear system now acts as shown by the superposition in Figure .21, not all the harmonics are so deeply attenuated that they become irrelavent as they were in the "high frequency" case. Thus, we get the output shown in Figure .22 which clearly is not sinusoidal. However, the fact that the harmonics are doubly attenuated (both by the coefficient effect and the linear filter effect) is seen by noting both that the output has lower mean value and is "more sinusoidal" than 217
PAGE 224
1 0.9 0.8 0.7 lfJ 0.6 ::> "' 0.5 ""' "' 0 0.4 "" 0.3 0.2 0.1 0 I I 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 FREQUENCY Figure 0.20: Frequency domain output of the Hammerstein nonlinear subsystem"low frequency" case with fundamental amplitude normalized. 218
PAGE 225
1 0.9 0.8 0.7 CFJ 0.6 t: z :::> 0: 0.5 "' "" 0 0.4 0.. 0.3 0.2 0.1 I I I I / / '. / \ : \ l \ l \ ' / \ ' I \ l \ / \ : \ \\. FREQUENCY Figure 0.21: Effect of linear systemfrequency domain. 219
PAGE 226
OUTPUT VS TIME FOR I!AMMERSTEIN SYSTEM "(\ 42 OJ t: z 40 "' .... 38 "' 0 36 34 v (\ (\ (\ (\ I v 0 50 100 150 200 250 300 TIME(SECS) Figure 0.22: Output of the Hammerstein system"low frequency" case. 220
PAGE 227
in the Wiener case. There is one further observation that should be discussed using this style of argument. Note that the difference between the mean output for the two models is greater for the high frequency case than for the low frequency case. This follows because the Taylor coefficients decrease with harmonic number, while the linear system attenuates with frequency. Thus, for the high frequency input case, the harmonics in the Hammerstein model are more completely removed by the linear filter effect than i n the low frequency input case, while the removal of harmonics by the decreasing coefficients of the Wiener model is the same whether the input frequency is low or high. This explains the greater difference in mean output for the high frequency case. 221
PAGE 228
BIBLIOGRAPHY [1] M. Aoki, State space modeling of time series, 2nd ed. SpringerVerlag, Berlin, 1990. [2] K. J. Astrom and B. Wittenmark, Adaptive Control. AddisonWesley, New York, 1989. [3] S. P. Banks, Mathematical theories of nonlinear systems. PrenticeHall In ternational, New York, 1988. [4] H. A. Barker, S. N. Obidegwu, and B. Pradisthayon, Performance of antisymmetric pseudorandom signals in the measurement of 2ndorder Volterra kernels by crosscorrelation, Proc. IEE. 119 (1972), 353362. [5] T. J. Barstow and P. A. Mole, Simulation of pulmonary 02 uptake during exercise transients in humans, J. Appl. Physiol. 63 (1987), 22532261. [6] D. A. Belsley, Conditioning diagnostics. John Wiley, New York, 1991. [7] J. S. Bendat and A. G. Piersol, Random data: analysis and measurement proceedures, 2nd ed. John Wiley, New York, 1986. [8] S. A. Billings and W. S. F. Voon, Leastsquares parameter estimation algorithms for nonlinear systems, Intern. J. System Science. 15 (1984), 601615.
PAGE 229
[9] S. A. Billings and W. S. F. Voon, A predictionerror and stepwiseregression estimation algorithm for nonlinear systems, Int. J. Control. 44 (1986), 803822. [10] G. A. Brooks and T. D. Fahey, Exercise physiology: human bioenergetics and its applications. John Wiley, New York, 1984. [11] J. Candy, Signal processing: the modern approach. McGrawHill, New York, 1988. [12] R. E. Crochiere and L. R. Rabiner, Multirate digital signal processing. PrenticeHall, Englewood Cliffs, N.J., 1983. [13] W. D. T. Davies, System identification for selfadaptive control. Wiley Interscience, London, 1970. [14] D. Essfeld and U. Hoffman, A model for studying the distortion of muscle oxygen uptake patterns by circulation parameters, Eur. J. A ppl. Physiol. 62 (1991), 8390. [15] G. F. Franklin, J.D. Powell, and M. L. Workman, Digital control of dynamic systems, znd ed. AddisonWesley, Reading, Ma. 1990. [16] Y. Fujihara, J. R. Hildebrandt, and J. Hildebrandt, Cardiorespiratory tran sients in exercising man, J. Appl. Physiol. 35 (1973), 5867. 223
PAGE 230
[17] S. A. Glantz and B. K. Slinker, Primer of applied regression and analysis of variance. McGrawHill, New York, 1990. [18] G. C. Goodwin and R. L. Payne, Dynamic system identification: expenmen tal design and data analysis. Academic Press, New York, 1977. [19] D. Graupe, Time series analysis, identification and adaptive filtering, 2nd ed. Robert E. Krieger Publishing, Malabar, Fl., 1989. [20] F. S. Grodens, J. Buell, and A. J. Bart, Mathematical analysis and digital simulation of the respiratory control system, J. Appl. Physiol. 22 (1967), 260278. [21] R. A. Haddad and T. W. Parsons, Digital signal processing: theory, appli cations, and hardware. W. H. Freeman, New York, 1991. [22] T. C. Hsia, System identification: leastsquares methods. D. C. Heath, Lexington, Mass., 1977. [23] R. L. Hughson, D. L. Sherrill, and G. D. Swanson, Kinetics of V02 with impulse and step exercise in humans, J. Appl. Physiol. 64 (1988), 451459. [24] N. L. Jones, Clinical exercise testing, 3rd ed. W. B. Saunders, Philadelphia, 1988. [25] R. H. Jones and N. C. Steele, Mathematics in communication theory. Ellis Horwood Limited, Chichester, England, 1989. 224
PAGE 231
[26] S.M. Kay, Modern spectral estimation theory and application. PrenticeHall, Englewood Cliffs, N.J., 1988. [27] I. J. Leontaritis and S. A. Billings, Inputoutput parameteric models for non linear systems, Int. J. Control. 41 (1985), 303344. [28] J. N. Little and L. Shure, The Signal Processing ToolboxTM user's guide. The Math Works, Inc., Natick, Ma., 1989. [29] L. Ljung, System identification: theory for the user. PrenticeHall, Engle wood Cliffs, N.J., 1987. [30] P. Z. Marmarelis and V. Z. Marmarelis, Analysis of physiological systems: the white noise approach. Plenum, New York, 1978. [31] W. D. McArdle, F. I. Katch, and V. L. Katch, Exercise physiology: energy, nutrition, and human performance, 2nd ed. Lea & Febiger, Philadelphia, 1986. [32] A. J. Miller, Subset selection m regresswn. Chapman and Hall, London, 1990. [33] Minitab, Inc., MIN IT ABfP reference manual. Minitab, Inc., State College, Pa., 1989. [34] R. R. Mohler, Bilinear control processes. Academic Press, New York, 1973. 225
PAGE 232
[35] M. B. Priestley, Nonlinear and stationary time senes analysis. Academic Press, San Diego, 1988. [36] W. J. Rugh, Nonlinear system theory: the Volterra/Wiener approach. Johns Hopkins University Press, Baltimore, 1981. [37] S. Sagara and Z. Zhao, Numerical integration approach to online identifica tion of continuoustime systems, Automatica. 26 (1990), 6374. [38] SAS Institute Inc., S AS/STAT user's guide, Version 6, Fourth Edition, Volume 2. SAS Institute Inc., Cary, N.C., 1989. [39] M. Schetzen, The Volterra and Wiener theories of nonlinear systems. John Wiley, New York, 1980. [40] N. B. Slonim and 1. H. Hamilton, Respiratory Physiology, 5th ed. C. V. Mosby, St. Louis, 1987. [41] T. Soderstrom and P. Stoica, System identification. PrenticeHall, Engle wood Cliffs, N.J., 1989. [42] S. D. Stearns and D. R. Hush, Digital signal analysis, 2nd ed. PrenticeHall, Englewood Cliffs, N.J., 1990. [43] H. J. A. F. Tulleken, Generalized binary noise testsignal concept for improved identificationexperiment design, Automatica. 26 (1990), 3749. 226
PAGE 233
[44] H. Undehauen and G. P. Rao, Identification of continuous systems. North Holland, Amsterdam, 1987. [45] B. Widrow and S.D. Stearns, Adaptive signal processing. PrenticeHall, En glewood Cliffs, N.J., 1985. [46] B. J. Whipp, S. A. Ward, N. Lamarra, J. A. Davis, and K. Wasserman, Parameters of ventilatory and gas exchange dynamics during exercise, J. Appl. Physiol. 52 (1982), 15061513. 227
