&EPA
United States
Environmental Protection
Agency
Environmental Research
Laboratory
Athens GA 30605
EPA 600 3 78-104
December 1978
Research and Development
New Sampling
Theory for
Measuring
Ecosystem Structure
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ECOLOGICAL RESEARCH series. This series
describes research on the effects of pollution on humans, plant and animal spe-
cies, and materials. Problems are assessed for their long- and short-term influ-
ences. Investigations include formation, transport, and pathway studies to deter-
mine the fate of pollutants and their effects. This work provides the technical basis
for setting standards to minimize undesirable changes in living organisms in the
aquatic, terrestrial, and atmospheric environments.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/3-78-104
December 1978
NEW SAMPLING THEORY FOR
MEASURING ECOSYSTEM STRUCTURE
by
Robert J. Mulhoi land
School of Electrical Engineering
Oklahoma State University
Stillwater, Oklahoma 74074
Grant No. R-805564010
Project Officer
James Hill, IV
Environmental Systems Branch
Environmental Research Laboratory
Athens, Georgia 30605
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GEORGIA 30605
-------
DISCLAIMER
This report has been reviewed by the Environmental Research Laboratory,
U.S. Environmental Protection Agency, Athens, Goergia, and approved for
publication. Approval does not signify that the contents necessarily reflect
the views and policies of the U.S. Environmental Protection Agency, nor does
mention of trade names or commercial products constitute endorsement or
recommendation for use.
-------
FOREWORD
Environmental protection efforts are increasingly directed towards
preventing adverse health and ecological effects associated with compounds
of natural or human origin. As part of this Laboratory's research on the
occurrence, movement, transformation, impact, and control of environmental
contaminants, the Environmental Systems Branch studies complexes of environ-
mental processes that affect the transport, transformation, degradation,
and impact of pollutants or other materials in soil and water.
This research includes the development and testing of mathematical
models that can predict general behavior of toxic organic materials in
terrestrial and aquatic ecosystems. The sampling theory and related systems
analyses in this report combine model development with preliminary testing
using small laboratory ecosystems. This method represents an explicit
procedure for designing laboratory ecosystem experiments that provide
appropriate tests of model validity.
David W. Duttweiler
Director
Environmental Research Laboratory
Athens, Georgia
-------
ABSTRACT
This research considered the application of systems analysis to the study
of laboratory ecosystems. The work concerned the development of a methodology
which was shown to be useful in the design of laboratory experiments, the pro-
cessing and interpretation of the results of these experiments, the develop-
ment of model structures for microcosms based upon experimental data, the
testing of the predictive capabilities of evaluative models, and the analysis
of measurement errors in the presence of noise. Development of the experi-
mental design methodology was based upon a non-Nyquist sampling theory which
gave a p^oii computable bounds on the sample period for data collection.
Results of the research project included a definition for the bandwidth
of an ecosystem in the frequency domain, an indication of the effect of mea-
surement noise on the sample period for data acquisition, a scheme for reduc-
ing this noise, and a computer code for model identification based upon the
experimental design developed. The project also produced the conditions under
which two new schemes for model identification were applicable to the con-
struction of linear compartment models. These new identification schemes were
based upon system input perturbations and state observations which were often
readily implemented in the laboratory.
The research objective was to produce a scheme for the evaluation of
microcosm models used for the assessment of the environmental transport,
transformation, degradation, and fate of toxic materials in aquatic systems.
An hypothesis was advanced that the evaluative scheme will be developed in the
near future into a procedure by which models are tested against laboratory
data and refined until they realistically represent perturbed natural eco-
systems receiving toxic inputs.
The report was submitted in fulfillment of Research Grant R-805564010 by
Oklahoma State University under the sponsorship of the U.S. Environmental Pro-
tection Agency. This report covers the period October 23, 1977, to October
24, 1978, and the work was completed as of November 1, 1978.
IV
-------
CONTENTS
Foreword iii
Abstract iv
Figures vi
Tables vii
Acknowledgments viii
1. Introduction 1
2. Conclusions and Recommendations 3
3. Sampling Theory and Applications .... 4
4. Frequency Domain Techniques 20
5. Identifiability 25
6. Constant Infusion 36
7. Experimental Design 45
8. Computer Program 51
References 55
-------
FIGURES
Number page
1 An example two compartment system with precursor 7
2 Probability density function for w. 11
3 Spectral graph of the estimated mean value for a,-, of Smith's
model '! 16
4 Spectral graph of the estimated variance for a,, of Smith's
model '! 17
5 Spectral graph of the estimated mean value for a,9 of Smith's
model '. 18
6 Spectral graph of the estimated variance for a,9 of Smith's
model 19
7 Effective bandwidth construction 21
8 Example Bode plot 23
9 Hypothetical compartment system 40
10 Diagnostic input for experimental design 46
11 Response to diagnostic input for experimental design 46
12 Flow chart for computer algorithm 52
13 Flow chart for subroutine EIVECT 53
VI
-------
TABLES
Number Page
1 Basic Microcosm Experiments 6
2 Truncated State Samples for T = 0.01 12
3 Chi-Squared Test for a,, with N = 400 and 17 Degrees of
Freedom X217;0>05 = 27.59 14
4 Chi-Squared Test for a,,, with N = 400 and 17 Degrees of
Freedom Y217 n nc = 27.59 15
l/;U.Ub
vn
-------
ACKNOWLEDGMENTS
It is a pleasure to acknowledge the constant support and sustained inter-
est in this research program provided by Dr. William R. Emanuel of the Envi-
ronmental Sciences Division of Oak Ridge National Laboratory. Specifically,
Dr. Emanuel provided the FORTRAN code for the subroutine EISPACK which proved
to be extremely useful for eigenvalue/eigenvector calculations.
It is often difficult to distinguish between basic and applied research.
And, with respect to the research described in this report, it is often hard
to discern the part funded by the EPA grant from that funded as basic research
in systems theory by NSF Grant ENG75-05341-A01.
viii
-------
SECTION 1
INTRODUCTION
This research report concerns an investigation into the application of
systems analysis to the study of laboratory ecosystems. Despite the recent
resurgence of interest in microcosms, there appears to be no clear understand-
ing of how to best use system theory and mathematical modeling in this impor-
tant research area. However, there is apparently some agreement as to the use
of laboratory ecosystems in ecological research. At least with respect to
general research objectives important to the U.S. Environmental Protection
Agency, the following reasons for microcosm research are mentioned by Draggan
(1976), Brockway, et al. (1978) and Van Voris, et al. (1978): (1) to develop
a screening method capable of assessing the environmental fate and effect of
chemicals, (2) to study basic ecosystem processes under controlled environmen-
tal conditions, (3) to study ecosystem level variables as measures for ecolog-
ical complexity, and (4) to test or verify mathematical models.
The appeal for microcosm studies is presumably based upon the replicabil-
ity of experiments which can be performed. This replicable nature of micro-
cosm experiments is generally assumed to result from the control which can be
exerted over the laboratory environment. Thus, an ecologist concerned with
the study of ecosystem function can employ the hypothesis that microcosms sub-
jected to similar physical and chemical conditions, even when inoculated with
samples derived from different sources, will develop into functionally similar
systems, see Brockway, et al. (1978). Two quantitative tools are required for
work based upon this hypothesis of functional invariance: (1) a rigorous ex-
perimental design for measuring ecosystem structure and (2) a theoretical con-
struct which relates ecosystem structure and function. The research reported
on here presents the idea that system theory and mathematical modeling can
provide both of these tools.
There are also diverse opinions on how to interpret or apply the results
of microcosm research. Many ecologists view microcosms as "stripped down" or
"simplified" biological models which are intended to simulate some specific
natural process (Odum, 1971, p. 21). On the other hand, Brockway, et al.
(1978) view laboratory ecosystems neither as biological models nor as analogs
of natural ecosystems. Instead, they take what may be a more realistic view-
point that microcosms are ecosystems in their own right which have their own
community structure, mineral cycling, and pathways for energy flow typical of
all natural systems. It can be expected that toxic material would disrupt
such integrated microcosm processes in a manner similar to any other ecosystem
of related structure and function. Thus, a mathematical model which is de-
signed to predict the fate and effect of a substance in any ecosystem which
represents a class of systems defined by specific functional attributes should
at least be capable of similar predictions for a microcosm belonging to the
same class. It follows that models developed for screening can be tested in
the laboratory against replicable results derived from perturbed microcosm
1
-------
experiments provided a rigorous theory for sampling ecosystem structure is
employed for data acquisition. Within the arithmetic precision of the experi-
mental data acquired from the microcosm, such a theory for sampling ecosystem
structure must completely capture the inherent dynamics of the mathematical
model. Thus, the use of system analysis is employed to help insure the
results of microcosm experiments are useful in testing the prediction of
ecological fate and effect made by mathematical models.
The idea that systems analysis has application to microcosm research is
not new. System applications usually focus on the use of models to explore
potentially interesting microcosm experiments, the testing of models against
microcosm data, and the use of basic model structure to define what variables
are of interest and hence what should be measured. The last potential use
mentioned may hint at the most important application for system analysis, i.e.
its use in the design of experiments. Mulholland and Gowdy (1977) presented,
for perhaps the first time, the concept that systems analysis and modeling can
be used to develop a quantitative theory for measuring ecological structure in
laboratory microcosm. The idea is to use systems analysis, in much the same
way statistics is now used, to design microcosm experiments. Starting with a
theory used to produce bounds on the temporal sampling period for measuring
compartment level variables in a microcosm, research during the past year has
produced a design for input function which places in evidence the complete
dynamic response resulting from a system level perturbation. The resultant
experimental design provides for the collection of sufficient data to identify
the parameters of two models types: (1) a sequence of (local) linear models
with a. pnJLonJL unconstrained structure and (2) a nonlinear model of postulated,
but not quantified, structure. The basic input function design also allows
for the collection of independent data useful in testing the responses of the
resultant models to perturbations.
In summary, research has disclosed an experimental design for laboratory
ecosystems which provides a method for introducing the appropriate system
level perturbations, introduced through natural system input ports, in order
to elicit full system dynamic response. The design also provides for ampli-
tude calibration and periodic time sampling which captures the full dynamic
state response. The state response data which result from the experimental
design establish a test for system linearity with respect to multiple pertur-
bations, a data base for constructing local linear models, an independent data
base for nonlinear model construction, and still another data base for testing
model response.
-------
SECTION 2
CONCLUSIONS AND RECOMMENDATIONS
The results of this research have produced an experimental design for the
laboratory analysis and evaluation of perturbed microcosms. The complete de-
sign has not yet been tested against actual experiments. A continuation of
this line of research would allow such a test to take place. Additional re-
search is also needed for the development of a stochastic version of the ex-
perimental design, thus combining statistics with systems analysis which could
lead to a powerful technique for error analysis. The development of a statis-
tical error analysis should give rise to techniques for filtering raw data,
and hence for minimizing measurement errors in the presence of noise. Testing
of mathematical models against microcosm derived data, based upon the proposed
development of a stochastic experimental design, could lead to a rigorous
methodology for model validation.
With respect to analysis in the frequency domain this research has been
unable to reconcile the Mulholland-Gowdy sampling theory with the Nyquist
theory used in engineering systems analysis. The main problem here, yet to be
resolved, concerns the development of an effective technique for bandlimiting
ecosystems. Unlike engineering systems, it may be impossible to define a use-
ful measure for the equivalent bandwidth of an ecosystem. The Mulholland-
Gowdy sampling theory provides bounds on the sample period used for data
collection in the laboratory, which in turn describe a window in the frequency
domain where the state response spectrum resides. However, it is unclear
whether or not this window, having been developed under prescribed input con-
ditions, gives a universal representation for the frequency spectrum of an
ecological system.
The study of two special forms of compartmental structure has produced a
generalization of the concept of identifiability. First, for systems, in
which the state is not directly observable, but measurements of only accumu-
lation pools are possible, the conditions for identifiability, controllability
and observability have been defined. Second, the conditions under which
systems receiving a constant infusion of labeled material are identifiable
have been produced. Although these conditions on the identifiability for the
constant infusion method are somewhat limited by system structure, it should
be noted that most system scientists believe that systems cannot be identified
from only steady-state data and that excitation of transient states are always
required. In spite of this misconception, it has been shown that complete
sets of model parameters can be obtained when systems are in steady-state
under specific conditions which are often met in practice. The study of
identifiability conditions is necessarily theoretical in nature. So, in order
to be of some use in experimental design, the conditions which have resulted
from this research must be tested in the laboratory against the data derived
from microcosm experiments.
-------
SECTION 3
SAMPLING THEORY AND APPLICATIONS
MODEL IDENTIFICATION
The determination of a model relating structure and function of a system
from experimental time series observations of the state is known as system
identification. It includes the estimation of system parameters for the ob-
servable components within the biological system. One of the fundamental
problems of system identification is whether or not the parameters can be
determined uniquely (Lee, 1964). If unique estimates cannot be obtained (the
system is not identifiable), then the experimental procedures should be al-
tered. Determination of structure from function results from an input-output
tracer experiment. The different identification schemes that are available
can be classified according to the basic elements of the problem and the input
signals, the earliest methods being those based on frequency, step and impulse
responses. Special inputs rather than those of normal operation are necessary
for these schemes, so the techniques are termed off-line identification and
therefore apply only to linear stationary processes where input/output rela-
tions which hold for one set of inputs are good for all inputs (Graupe, 1972).
For the purposes of tracer analysis the local model is a linear donor
controlled biological compartment model defined by
x^ = Ax + Eto (3.1)
where the state x_ is an n-dimensional vector of the compartmental contents,
the coefficient matrix A_ = (a..) is an ordered array of the compartmental rate
coefficients, the input ^defines the ecosystem inputs as an m-dimensional
vector, and the input connection matrix J3 = (b..) is of n rows and m columns
with entry b. . equal to one if input u- enters compartment i, if not, the
I J J
entry is zero.
In equation (3.1), the inputs u^ and their entry into the ecosystem com-
partments, the B^matrix, are known. The solution of the identification prob-
lem for (3.1) requires a method to determine the entries in A, containing the
rate coefficients, given sufficient measurements of the state x_.
BASIC EXPERIMENT
The use of radioactive tracers is prone to varying experimental proce-
dures, especially with the advent of highly sophisticated instruments and new
radioisotopes. Practical experiments must be designed to coexist with devel-
oped theoretical techniques. A standard method has developed using three
basic experiments to assess the environmental impact on aquatic microcosms of
-------
toxic material inputs. As shown in Table 1, the first test input represents a
slug of toxic material, the second considers a constant chronic input level of
toxin, and the third forms the superposition of the aforementioned experi-
ments. The microcosm responses to the three basic inputs are also shown in
Table 1. The toxic slug of part (a) elicits a transient response with initial
condition XQ, and the constant input level of part (b) gives rise to a steady-
state defined by x . By assuming a linear model for the microcosm system,
part (c) of Table 1 illustrates the resultant superposition principle. It
should be noted that the basic microcosm experiments contain information re-
garding the validity of the linear modeling assumption. For linear models,
the microcosm responses to the inputs of Table 1, parts (a) and (b), should
sum to give the response to the input of part (c). If this superposition is
in fact true, the single dose method [the input of part (a)] can be effective-
ly used to identify the numerical model for the microcosm. Therefore, the
validity of the linear state model (3.1) should be clear from the outcomes of
the basic experiments of Table 1. Generally, a local linear model can be
identified for small perturbations (x is small compared to x ) caused by im-
pulse inputs at a given system loading defined by u . A local model defined
in this way will be dependent upon the system loading.
Also in Table 1, two new experiments are proposed for linear model iden-
tification using tracers or toxic material flows as microcosm inputs. These
inputs, based upon step functions, have been designed for ease of application
in the laboratory and mathematical convenience for system identification.
These basic experimental inputs lead to a sampling theory which is essential
to the successful application of standard identification schemes. This re-
sult, in turn, prescribes an experimental design for the a p^ioii determin-
ation of the period for data collection in the laboratory.
The basic experiment of part (e) of Table 1 has been designed as an iden-
tification technique for local linear models. Linearity can be tested in the
laboratory by observing response superposition as result of adding two micro-
cosm inputs.
As previously noted, the basic experiments (d) and (e) of Table 1 have
been designed for ease of application in the laboratory and to aid in the de-
velopment of a sampling theory. It turns out that under certain conditions
these basic experiments also avoid the problem of tracer analysis of having
input access to each compartment. However, the equivalence between the inputs
of parts (c) and (e) of Table 1 is essential in eliminating the so-called pre-
cursor problem. In Figure 1 an example of a two compartment model with a
single dose input into the precursor compartment is shown in part (a) along
with the input and compartment response plots in part (b).
The precursor problem deals with the impossibility of estimating the max-
imum value (x ) of the response of compartment x? to the impulse input to com-
partment x,. In order to develop an effective laboratory data collection
strategy, it is necessary to have known rules for amplitude and time discret-
ization. Because x is known before the experiments is completed, it is very
difficult in ordinary tracer analysis to develop a rule for amplitude quanti-
zation. Given a sample period, data on microcosm structure is collected at
discrete times and the values of these data are quantized by the finite arith-
metic precision of the measuring instruments. However, the quantization of
-------
TABLE 1. BASIC MICROCOSM EXPERIMENTS
Input
Response
(a)
(b)
(c)
(XS
XS
U
(d)
t \ U?
(e) Z
-*
t
X<
X,
-------
1 output
impulse
input
X1
transfer _
«,
output
Figure 1. An example two compartment system with precursor.
time and amplitude are not sequential processes; they are, in fact, related to
one another. This gives rise to the commonly used laboratory procedure for
single dose experiments of taking microcosm samples very rapidly at first
until Xr, reaches its maximum and then more slowly to accommodate the compart-
mental dynamics inherent in the free response. The asynchronous sampling
procedure described is difficult to analyze, particularly when the propagation
of measurement errors through the identification algorithm is required. A
periodic sampling scheme would be more desirable.
The basic experiments (d) and (e) of Table 1 eliminate the precursor
problem. This follows from the fact that transitions between steady-states in
compartmental system are monotonic and bounded by the initial and final states
(Thron, 1972). Therefore, the amplitude quantization process is well defined,
and a periodic sampling strategy is possible.
SAMPLING THEORY
The collection of time series data for each compartment x.(t) in the bio-
logical system takes place in two distinct steps in which first time and then
amplitude are quantized. The samples of the state are obtained only at dis-
crete times. These samples are periodically spaced apart by T time units
giving a data base defined by 1x^(0), X(T), ... , x[(k+l)x]}. It has been
shown that the minimum data base for solution of the identification problem is
(n+1) state samples or k = n. Next, the finite precision with which the
measurements of state are made must be accounted for in terms of a mathemat-
ical relationship. The assumption is that x.(t) and its samples are fixed
point numbers which are less than or equal to one. These samples, arising
from radiological measurement equipment, are data with no more than a fixed
number of digits to the right of the decimal point. Let 10" , where d is a
positive integer, be the smallest non-zero number in the data set.
-------
Because the eigenvalues of the system described by (3.1) are bounded by
Gerschgorin circles of known radii (Hearon, 1963), it is not difficult to show
that transitions between steady-states, as in Table 1, are bounded by two ex-
ponential factors (see Mulholland and Gowdy, 1977):
exp(-at) < x.j(t) <_ exp(-fit) (3.2)
where a/2 = T is the maximum turnover and 6 is the minimum non-zero trans-
fer out of the system. The selection of the sample period T for data collec-
tion is now possible. Using (3.2), it is easily shown that T is bounded by
T~ <_ T <_ T , where T is the maximum sample period and r~ is the minimum,
both of which are determined by the arithmetic precision of the state samples.
It is possible for the sample period to be too small. This occurs when
x.j(t) and x^t + r) have the same arithmetic precision; that is, when x.(t + t)
is truncated to d digits, it equals the truncated value of x.(t). Since (3.2)
gives an upper bound for x^(t), a measure of the minimum sampling time, de-
noted by T~, is given by
l-10"d = exp(-6O < XI(T") (3.3)
or solving for the minimum sampling period yields
T = (-1/6) In (l-10"d). (3.4)
If T > T~ the state samples are independent, but an upper bound on T
exists. Suppose T is chosen large enough so that x. (rvr), the last sample in
the data set, is truncated to zero, then the data set is incomplete and the
identification problem cannot be solved. Thus, the maximum sample period,
denoted by T , must be such that the (n+l)-th sample is at least equal to the
smallest number in the finite precision arithmetic; that, x.(nr ) 21 10~ But
(3.2) gives a lower bound for x.(t); thus, the maximum sample period is easily
computed from
10"d = exp(-anT+) < x^m*), (3.5)
which gives
T+ = (-l/ZnT^1) In (10"d). (3.6)
Based upon intuition, experience, isolated laboratory experiments, and/or
estimates from the literature, biologists tend to know the turnovers for the
compartments present in a system. Such knowledge is at least accurate to an
order of magnitude. Thus, the maximum sample period is easily determined from
(3.6).
8
-------
IDENTIFICATION ALGORITHM
Consider the experiment of part (d) of Table 1 in which a constant
(known) input is applied to an ecosystem modeled by (3.1). After sufficient
time has elapsed, a steady-state is observed, as indicated by x^ = 0_. Denote
this steady-state by the vector x , which is the solution of
0 = Ax.
Bu.
(3.7)
In the experiment under discussion, the input is assumed to be a constant rate
of infusion of labeled material. It is also assumed the steady-state can be
measured using standard radiological methods. On a new time scale, the exper-
iment is re-started at time t = 0 by removing all inputs to the ecosystem.
Thus, before t = 0, the system is in steady-state x_(0) = x^, while for t > 0
the state of the system is governed by the solution of
= Ax
(3.8)
because u_ = 0. for t > 0. The labeled contents of the compartments of the
ecosystem start at t = 0 with values determined by x^ and decrease exponen-
tially as t becomes large as all the labeled material flushes from the system
along the pathways for material transfer.
Discretized in time the solution of (3.8) is governed by
x (IT) = * X[(I-I)T] (i = 1, 2, ... ),
(3.9)
where $_ = exp(Ar) with T~ £ T £ T .
A continuous record of the state x_(t) is not required in order to deter-
mine $_, which by definition gives A_. Since £ is an n x n constant matrix, con-
2
straint equations numbering n are required to prescribe the entires in this
matrix. The problem reduces to finding A from a finite number of discrete and
regularly spaced samples of the state, given at times t = 0, T, 2x, ... , kr
where T is the sample period. Substitution of these state samples into (3.9)
gives
X(T) = $.x_
x.(kr) = ix_
(3.10)
j
This linear algebraic system contains n unknowns (all the entries of £) and k
independent vector equations of n components each, fora total of kn equations.
For k < n no unique solution of (3.10) exists, and the system is overdeter-
mined for k > n. Thus k = n or (n+1) discrete samples of the state are neces-
sary and sufficient to determine £.
-------
It is possible to rewrite system (3.10) in the following more standard
algebraic form
X2 = *Xr (3.11)
where X^ is an n x k matrix formed by columns from the first k state samples,
, X_(T), ... , x_[(k-l)T]>, (3.12)
and X_2 is similarly obtained from the next k state samples delayed by one
sample period,
), x_(2x), ... , x.(kT)}. (3.13)
For k = n in (3.10) the solution for £ is obtained by inverting the matrix X,.
This assumes n linearly independent state samples which is indeed the case
when the ecosystem model defined by (3.1) applies. Thus, the identification
problem is solved by
A = (1/T) ln(X_2X^ ') (3.14)
under these stated conditions.
It is mathematically necessary and sufficient to make (n + l) discrete
state measurements in order to solve for $_. At least (n+l) measurements are
required for a unique solution, however on occasion it may be desirable to
obtain more than (n+l) state samples. For example, when significant experi-
mental error is suspected in the state data, it may be of use to increase the
dimension of the data set beyond (n+l). Under these conditions k > n, and
there are more equations than unknowns in (3.10). A Gaussian least squares
estimate for the solution of A can then be obtained (Lee, 1964).
SIMULATED EXPERIMENT
Early in the project, it was surmized that problems might arise in get-
ting the data from the EPA planned microcosm experiments which were designed
to test the Mulholland-Gowdy sampling theory. So, a simulation of experimen-
tal data was devised using an ecosystem model. It is difficult to theoretic-
ally analyze the propagation of data noise through a nonlinear algorithm.
Since the identification algorithm is nonlinear, a simulation was used to
examine the effects of data noise errors (on the algorithm). In the simula-
tions, data noise was propagated through the algorithm to obtain the statis-
tics of individual parameters within the identified model.
The identification algorithm,defined by (3.11), requires discrete samples
of the system state: x^ = x^kx). Actually, the ensemble of these samples is
itself a sample function for a stochastic process. Because of noise inherent
in the microcosm processes and errors introduced by the measurement techniques
and instruments, the analysis of the basic experiment (defined in Table 1)
cannot be wholly based upon deterministic systems theory, but it must include
a statistical analysis. The statistical interpretation of the basic microcosm
experiment relies upon the fact that the initial and final steady-states are
random variables with measurable statistics. These steady-states are connec-
ted by a sample function of a stochastic process which is a statistical repre-
10
-------
sentation of the sequence of transition states.
The stochastic process which represents the transition states has been
simulated by a deterministic variable (z, ) perturbed by a uniformly distrib-
uted random variable (w, ).
\ = \ + ^k (3'15)
for each k = 0, 1, ... , n. The motivation for this representation derives
from knowing (by measurement) the system state to an arithmetic precision of
10" . Then, the random variable w. is uniformly distributed from 0 to 10
with a probability density function given by Figure 2. This perturbation
simulates noise superimposed upon the truncated state.
io-d
Figure 2. Probability density function for w. .
The solution of (3.11) for the identified model shows in (3.14) a depen-
dence upon the sample period. Indeed, there probably exists an optimal sample
period which depends upon the statistics of the experiment. In order to test
for this optimal sample period and in general study how noise propagates
through the identification algorithm to the entries in the model A-matrix, a
Monte Carlo experiment was decided upon. The exact solution of the model,
correct to d-digits, has been obtained. The state is then truncated to d-
digits and corrupted by noise according to (3.15). Simulated noise measure-
ment data, defined by (3.15) are passed through the identification algorithm
many times until the statistics of the entires in the A-matrix are measure-
able. The following example will clarify this scheme.
The ecosystem model chosen to represent, through computer simulations,
the effect of measurement noise on the identification algorithm and sampling
theory is a hypothetical ecosystem first described by Smith (1970). Smith's
system is not a model of a real system, but is a model built upon realistic
numerical parameters and a state response which mimics an aquatic ecosystem
operating under a phosphorus loading. Very simply, Smith's system consists of
three compartments: water (x,), an aquatic plant population (x2), and an
herbivore population (x-J. The only ecosystem input (u,) is the phosphorus
contained in the water flowing through the system. The outputs are phosphorus
from the water compartment and from the herbivore compartment due to emigra-
tion or emergence. The plants exchange phosphorus with the water compartment,
11
-------
which is transferred to the herbivores through grazing.
From data presented by Smith, the rate coefficients can be computed,
which results in
x-j = x2 + x3 - x1 + u,
x2 = 14x1 - 95x2 (3,16)
x3 = 90x2 - 14x3
as a (local) linear donor controlled model defined near the prescribed input
of u.j =100 mgP/day. The time scale of the model is defined by the rate co-
efficients which are given in terms of units of reciprocal days.
With u1 = 100, equation (3.16) was simulated on the IBM-370 using the
CSMP simulation language package. State samples correct to 4 1/2 decimal
places, i.e. x2 takes initial values with 5 significant digits but it eventu-
ally declines to 4 signigicant digits, were obtained for 20 different sample
periods ranging between 0.005 and 0.120. It should be noted that (3.16) with
a = 190 and 6 = 2 has sample period bounds of T~ = 5 x 10 and T+ = 0.0202
for d = 5, and T~ = 5 x 10~ and T+ = 0.016 for d = 4. For T = 0.01, the four
state samples are given in Table 2. These samples illustrate the particular
basic experiment performed where the initial steady-state entrained by the
constant input is the initial condition. The input is then removed and the
state declines, according to the system dynamics, monotonically to zero.
TABLE 2. TRUNCATED STATE SAMPLES FOR T = 0.01
k
0
1
2
3
9
8
7
7
xl
.5000
.5748
.7807
.0939
1
1
1
1
x2
.4000
.3506
.2549
.1519
9
8
8
8
X3
.0000
.9844
.9052
.7513
A pseudo-random number generator was then used to drive the Monte Carlo
simulations for model identification. The density function for w. , uniformly
-4 -5
distributed from 0 to 10 , gives a mean y = 5 x 10 and a variance of
j -in w
a = 8.33 x 10 . Using standard statistical tests with a level of signif-
W
icanceof 0.05, it was determined that 400 samples were required to estimate
12
-------
the mean and variance for w. . It was then assumed that 400 computer runs of
the identification algorithm were sufficient.
For each sample period, the random number w. was added to the simulated
states z. to form the measured data x,, for the identification algorithm. A
Ix (x
total of 400 models were generated for each T. Using the estimates given by
N
ya = (1/N) E a. (3.17)
a i=l 1
and
-9 N 9
a/ = (1/N-l) E (a, - y_r (3.18)
a i=1 i a
for N = 400, the mean and variance for the a,, and a,,, entries in the A-matrix
of (3.16) were determined.
A chi-squared fitness test was performed to test for normality of the
distributions for a,, and a-.^- The results are given in Tables 3 and 4.
Spectral plots of y and 3 for a,, and a,0 appear in Figures 3 through
a a ii i c.
6. For this particular model ecosystem, the empirical data suggest bounds of
0.007 and 0.07 for the sample period T as compared to T~ = 5 x 10 and T =
0.02 as computed from the theory for d = 5. Bounds for T cannot be specified
more closely because of the statistical nature of the data. For example, the
estimated mean for a,, is always within acceptable limits of the true value of
a-|, = -16 (see Figure 3) regardless of the sample period and independent of
the estimated variance. Therefore, the error inherent within the identifica-
tion algorithm cannot be exactly specified. However, the spectral plots in-
dicate a relationship exists between the sample variance and the sample per-
iod. It is possible that future research efforts will disclose a statistical
theory for this relationship.
Finally, a comment on the apparent disagreement between the theoretical
sample period bounds and those derived empirically from this specific example
seems to be in order. The theoretical bounds given in (3.4) and (3.6) are
general bounds for all compartmental systems with maximum turnover of a/2 = 95
and minimum non-zero exogenous turnover of 6 = 2, and that these bounds are
exceeded for a specific model does not invalidate the theory. Indeed, a
sample period chosen less than T = 0.02, say T = 0.01, appears to be perfect-
ly acceptable for model identification.
13
-------
TABLE
N = 400 AND
3. CHI-SQUARED TEST FOR a^
17 DEGREES OF FREEDOM X2-.7.0
WITH
.05 = 27'59
T
.005
.006
.007
.008
.009
.010
.011
.012
.013
.020
.030
.040
.050
.060
.070
.080
.090
.100
.110
.120
5a
-16.024
-15.973
-15.964
-16.057
-16.039
-16.010
-16.010
-15.983
-16.022
-16.016
-15.992
-16.050
-15.997
-15.965
-16.007
-15.780
-16.249
-16.190
-16.178
-16.368
o
aa
1.39 x 10"2
7.90 x 10"3
4.98 x 10"3
3.43 x 10"3
2.51 x 10~3
1.95 x 10"3
1.57 x 10"3
1.31 x 10"3
1.12 x 10"3
6.38 x 10"4
6.39 x 10~4
9.89 x 10"4
1.86 x 10"3
3.95 x 10"3
9.67 x 10~3
2.02 x 10"2
2,14 x 10"2
3.40 x 10"2
5.57 x 10"2
9.33 x 10"2
a(chi )
32.4
31.4
36.2
35.3
41.7
41.2
31.1
25.7
23.6
8.2
7.0
13.7
25.2
10.8
15.1
20.9
6.8
14.5
19.7
55.7
14
-------
TABLE 4. CHI-SQUARED TEST FOR a]2 WITH
N = 400 AND 17 DEGREES OF FREEDOM X^y.g Q5 = 27'59
T
.005
.006
.007
.008
.009
.010
.011
.012
.013
.020
.030
.040
.050
.060
.070
.080
.090
.100
.110
.120
/N
^a
5.147
4.793
4.541
5.578
5.404
5.035
5.184
4.888
5.165
5.113
4.947
5.322
4.984
4.786
5.046
3.680
6.516
6.160
6.081
7.306
- 2
aa
2.655
1.194
0.627
0.372
0.241
0.168
0.124
0.096
0.077
0.033
1.028
0.040
0.076
0.145
0.345
0.701
0.781
1.262
2.096
3.485
a(chi )
27.2
23.3
20.6
19.2
21.2
18.2
16.5
16.9
13.9
5.7
20.8
14.3
22.8
14.4
12.5
19.5
9.1
12.9
17.7
45.1
15
-------
-14.0
-15.0
-16.0
en
-170
-18.0
JRUE
VALUE
,1
.001 .003 .01 .03
.3
Figure 3. Spectral graph of the estimated mean value for a,, of Smith's model.
-------
.10
.08
.06
.04
.02
0
.001
.003
o
.01
.03
.3
T
Figure 4. Spectral graph of the estimated variance for a,, of Smith's model.
-------
7.0
6.0
00
4.0
I i I i i i r
I r
3.0
.001 .003 .01 .03
TRUE
VALUE
.3
Figure 5. Spectral graph of the estimated mean value for a,2 of Smith's model
-------
3.0
2.0
1.0
.8
J I
01 .003 .01 .3
r
.3
Figure 6. Spectral graph of the estimated variance for a,2 of Smith's model
-------
SECTION 4
FREQUENCY DOMAIN TECHNIQUES
For both ecosystem modeling and experimental data analysis, the tech-
niques of frequency domain analysis, using the Fourier series and transform,
are becoming more popular with ecologists, e.g., see Platt and Denman (1975),
Waide et al. (1974), Shugart (1978), Emanuel et al. (1978), and Van Voris et
al. (1978). This is intriguing because the common method of analysis for
determining the sample period for a communication system involves the frequen-
cy domain concept of bandwidth, expressed in terms of the Nyquist theory.
Therefore, it would be useful to have a system bandwidth defined for the sam-
pling and quantizing precesses used in the solution of the identification
problem.
A transfer function can be defined for linear systems with one input u(t)
and one output x(t). When the single input equals the unit impulse, also
called the Dirac delta function 6(t), the output is called the impulse re-
sponse and is denoted by h(t). The transfer function for the system is then
defined as the Fourier integral transform of the impulse response:
H(f) = _mr h(t) e"1a)tdt (4.1)
where H(f) is the transfer function, f = w/2ir is the frequency independent
variable measured in reciprocal time units and i = AT.
A sequence of n transfer functions for an ecosystem can be defined by
choosing, in turn, each compartment x. (j = 1, ... , n) of the state vector x^
as an output. The n transfer functions are then calculated by using
Hj(f) = Qr Xj(t)e-ia)tdt (4.2)
for all j = 1, ... , n.
Transfer functions reflect the dynamics of a system through the concept
of bandwidth, loosely defined as the range of frequencies which the system
described by H(f) passes in a relatively unchanged form. Carlson (1975) notes
that the rate at which a system can vary stored energy is determined by the
frequency response to perturbations as measured by the bandwidth. Thus, for
ecosystems, bandwidth and compartmental turnover appear to be related.
The so-called bandlimited systems have transfer functions which satisfy
H(f) = 0, for all f > B. (4.3)
For such systems the bandwidth is easily defined to be equal to B. However,
most systems, including ecosystems, are not bandlimited, so other methods of
defining the bandwidth are employed. One such method, called by Carlson
20
-------
(1975) the effective bandwidth, uses the equal area concept. Consider Figure
7 which plots and compares the amplitude of a typical transfer function (the
solid curve) with that of a rectangular window (the dashed curve). The rec-
tangular window is clearly bandlimited with bandwidth of B. The effective
bandwidth for H(f) is equal to B under the constraint:
2B | H(0) | = _af | H(f) | df
which amounts to finding B to give equal area under the two curves.
H(f)|
(4.4)
ja i
t
-BOB
Figure 7. Effective bandwidth construction.
It is now possible to compute the bandwidth for the ecosystem transfer
functions defined by (4.2). First, recall that x.(t) is amplitude limited by
(3.2). Next, let f = 0 in (4.2), giving J
Hj(0) =
(4.5)
Substituting (3.2) into (4.5), gives the following inequality
I/a < H.(0) < 1/6 (4.6)
J
which holds for all j = 1, ... , n. In order to complete the computation of B
from (4.4), it is necessary to evaluate the integral which appears in this
equation. Estimating this integral is possible through the definition of the
inverse Fourier transform:
f. (4.7)
(4.8)
XjU) - _00/°°Hj
Taking absolute values of (4.7) for all t > 0 and j = 1, ... , n, yields
H,(f)
\J
| df.
It is not possible to fully combine the inequalities (4.6) and (4.7) into
(4.4) to estimate the bandwidth, the best result being 1 <_ 2B/6. Thus, the
bandwidth of an ecosystem depends upon turnover, but not in the ordinary
sense. This result implies the bandwidth B to be greater than or equal to
one-half the minimum exogenous turnover, that is the transfer rate of material
out of the system to the environment.
For bandlimited frequency domain representations, the Nyquist sampling
theorem gives T £ B/2, where T is the sample period from which measured data
21
-------
can be used to reconstruct continuous time functions of system responses. The
two inequalities derived are incommensurate, however assuming 6 = 2B as a rule
of thumb, the Nyquist result implies that ecosystems should be sampled with a
sample period less than or equal to the reciprocal of the minimum exogenous
turnover. Another rule of thumb is obtained by assuming equality in (4.8) and
H.(0)a = 1 in (4.6) giving 2B/a = 1. Since a/2 = T -1, the bandwidth of an
J 0
ecosystem is then equal to the maximum compartmental turnover.
For systems which are not bandlinvited the bandwidth concept is at least
elusive. The bandwidth of an ecosystem has not been clearly defined, there-
fore it may be useful to consider the results presented in the light of an
example. As an elementary example of the frequency domain techniques pre-
sented, consider again Smith's hypothetical ecosystem defined by (3.16).
The transfer function H.(f) defines the frequency response of compartment
J
j with respect to the input which enters compartment one. The Bode plot for
the amplitude response is obtained from
Logarithmic gain = 20 log,JH.(f) |, (4.9)
< v J
which is measured in units of decibels (db). Plots of (4.9) for the three
compartmental transfer functions of Smith's model are given in Figure 8. These
Bode plots of the system transfer functions, each of which contains several
poles and zeros, are obtained by adding the individual plots due to each pole
and zero. Furthermore, this procedure is simplified by using only the asymp-
totic approximations to these curves, as shown in Figure 8. These Bode plots
describe the frequency response of the example ecosystem.
From (3.16) simple calculations yield 6=2 and a = 190. Therefore,
using T = /T"*"T~ and (3.6) and (3.4), the sample period is T = 0.00078 day for
d = 4, which can be rounded to T = 0.001 day.
As for the ecosystem bandwidth, the inequality 1 £ 2B/6 gives B >^
1(I/day). Figure 8 indicates the bandwidth to be approximately equal to 100
(I/day). Equating T = T, results in B £ 500 (I/day). Thus, the inequality
1 <_ 2B/6 yields very little useful information regarding the bandwidth of this
ecosystem, but the rule of thumb, derived as 2B/a = 1, gives B = 95 (I/day)
which seems reasonable.
Consider Figure 8 and measure the logarithmic gain at f = 95 (I/day) for
each compartmental transfer function. Note that |H.j(95)| > |H2(95)| >
|H3(95)|, which implies B-j > B2 > B3 as an ordering of the bandwidths for the
transfer functions of the three compartments. This ordering in turn implies
T, < T0 < T- for the three compartmental turnover times where T. = I/a...
1 L. J J J J
The ordering of the turnover times given is ecologically reasonable. In this
open (flow through) system, the phosphorus dynamics of the water compartment
are fast, more rapid than the uptake and transfer of phosphorus by the plants.
Also, the dynamics of phosphorus in the plants should be faster than for ani-
mals, as indicated. All this tends to add ecological intuition to the inter-
pretation of the mathematics of determining the ecosystem bandwidth using B =
V1-
22
-------
Odb
-60db-
1000
Figure 8. Example Bode plot.
23
-------
For systems which are not bandlimited the frequency domain concept of
bandwidth is inexact. Thus, the bandwidth of an ecosystem is not clearly
defined. However, based upon the calculations presented in this paper, B =
TQ , taken as an approximation, appears to be useful in expressing the eco-
system bandwidth as the maximum compartmental turnover.
The Mulholland-Gowdy theory for the periodic sampling of ecosystem
structure gives
T" < T < T , (4.10)
where the upper (T ) and lower (T") bounds on the sample period T depend on
the system turnovers, the arithmetic precision of the data, and the number of
system compartments. Processing of the transition state data, which are of
bandpass character, requires a window in the frequency domain given by
f& 1 f £ fu, (4.11)
where f = I/T~ and f = I/T describe, respectively, the highest and lowest
u x>
frequency (f) present in the data. If an assumption is made that the deter-
ministic system is perturbed by a white noise Gaussian random process, which
contains all frequencies with equally probable magnitudes, then the high fre-
quency noise can be filtered from the experimental measurements by bandlinvit-
ing data plus noise to an upper frequency bound given by f = I/T". Although
the transition state data is bandpass, the steady-state data is low pass re-
quiring measurements at essentially zero frequency, hence the low frequency
noise cannot be filtered from the data as easily as the high frequencies.
Least squares methods of data analysis remain the best solution to minimizing
low frequency noise. It is important to study the effect of filtering on the
error analysis of the data plus noise which come from the microcosm experi-
ments. _ _4
For Smith's hypothetical ecosystem the value for T" is 5 x 10" assuming
6 = 2 and d = 4. This gives a value of fy = 2000 per day, for which the Bode
plots for all three compartments are well below -60 db. Thus, a low pass
filter with an upper cutoff frequency of 2000 per day would surely pass the
information content of the microcosm experiments while rejecting the high
frequency noise present in the data.
24
-------
SECTION 5
IDENTIFIABILITY
Before implementation of the identification algorithm, conditions should
be investigated that will insure the uniqueness of a numerical estimation of
the associated parameters. A ptu.o/u. knowledge of which and/or how many param-
eters of the flow matrix A_ can be determined by means of the chosen input-out-
put experiment is termed identifiability (Cobelli, 1976). If the internal
couplings cannot be uniquely determined, changes in the mathematical model
structure or the experiment itself must be modified. Therefore, an identifi-
ability analysis may be considered as a necesseary step preceding identifica-
tion.
Often, identification is derived from numerical values of a few known
parameters and is reduced to a problem of state estimation in which it is
natural to deal with identifiability in a probabilistic manner. Identifiabil-
ity in such stochastic systems may be defined by the existence of an estimator
(Tse, 1973, 1978) or by the uniqueness of the probability distribution func-
tion (Bowden, 1973). A survey of similar investigations is provided by Reid
(1977). In the study of biological systems, identification often depends upon
a pfviofu, knowledge of the model structure (at least the order) and it is
appropriate to use a deterministic canonical system approach to the determin-
ation of model parameters. Other approaches to the problem include paramet-
rization (Glover, 1974) and parameter sensitivity analysis (Reid, 1977), which
are both closely related to transfer function methods.
An explicit discussion of the problem of structural identfiability in the
frequency domain is given by Bellman and Astrom (1970) where they propose the
use of the transfer function to investigate identifiability. Consider the
system
x(t) = A x(t) + B u_(t) (5.la)
y_(t) = C_,x(t) (5.1b)
which is the linear compartment model defined by (3.1) with a single dose
input u[(t), where y_(t) is an m-vector of measurements of the state x^ and C_ is
a p x n connectivity matrix of the form
C.- = o (5.2)
V
and
25
-------
for all i = 1, ... , p. Let Y_(s) be the Laplace transform of y(t) and U(s)
that of u_(t). Then ~
Y(s) = G(s) U(s) (5.3)
where Gjs) is the transfer matrix given by
G.(s) = C_(sl_ _ A)"1 B. (5.4)
and l_ is the n-th order identity matrix. The problem of identifying the flow
matrix A_ from the experimental data may be reduced to that of identifying the
coefficients of the transfer matrix, each coefficient being a combination of
the elements of A. Therefore, the A matrix is identifiable when this non-
linear system of equations has a unique solution.
Lin and Yu (1977) considered the single input and single output special
case of (5.1) in canonical form, where A is the companion matrix associated
with the characteristic polynomial
sm + am_1 s"1'1 + ... + alS + aQ = 0 (5.5)
and
^ = [c]c2 ... cn] (5.6)
B1 = [00... 1]. (5.7)
This leads to a scalar input/outout formulation for (5.3) in which G(s) con-
tains 2n unknown parameters (n poles, n-1 zeros, and a gain factor). The
identification of these parameters is demonstrated by a time series analysis
of the output variable given by
xlt V
y(t) = YI e + ... yn e .
A numerical algorithm called Prony's method (Hildebrand, 1974) is applied to
y(t) in order to resolve the parameters X. and y,- of the exponential curve
fit. The X. and y,- parameters are then shown to uniquely determine the 2n
unknown parameters of G_(s).
The precise definition of identifiability provided by Bellman and Astrom
(1970) is equivalent to saying that the parameters of A are identifiable if
and only if any change in the values of the coefficients causes a change in
the measured quantities y_(t). Thus, a necessary condition for identifiability
is that a system be controllable and observable in the sense of Kalman (1963).
Definition 5.1
A linear system is said to be completely state controllable if for any
initial state at t it can be driven to any other state in the state space at
t = T. °
26
-------
Definition 5.2
A linear system is said to be completely output controllable if it is
possible to drive any initial output at t to any final output in a finite
time interval.
Definition 5.3
A linear system is said to be completely observable if for every initial
time t and some finite time T,, t < t < T, every initial state at t can be
u 10 i o
determined from the knowledge of the output on t < t < T.
0
Rigorous treatments on the properties of controllability and observabil-
ity can be found in system theory literature (Kalman, 1963; Luenberger, 1963;
Lee, 1964; Ogata, 1967; and Chen, 1970). In the following, these properties
will be considered as they apply to the previously discussed methods of
identification for compartmental models.
DIRECT METHOD
Assume that a single dose of tracer may be directly inserted via special
routes into system compartments and that tracer behavior can be experimentally
monitored in specified compartments when the system is switched to a tracer-
free environment. The controllability and observability of such a system has
been explicitly discussed by Bellman and Astrom (1970) and Johnson (1976).
A system described by (5.1) is controllable if and only if the column
vectors of the n x nm test matrix
P = [B; AB; A2B; ...; A11"1 B] (5.8)
__£
span the n-dimensional space on which x^(t) is defined. In other words, P
must be of rank n if the system is to be controllable (Kalman, 1963).
Definition 5.4
A system is open (closed) when there is some (no) exchange of material
from the system either to the environment or to another subsystem. Using the
properties of the controllability matrix P. Johnson (1976) shows that closed
systems are uncontrollable regardless of their structure. Open systems are
also examined and are found to be controllable, in general, if the flow matrix
A is of rank n or (n-1) and at least one flow between the system and its
environment is controlled. Necessary conditions for complete controllability
can also be derived from physical reasoning. Cobelli and Romanin-Jacur (1976)
adopt a structural point of view by relating controllability and observability
of compartmental systems to their experimental design. A system is shown to
be completely controllable only if there exists at least one path from some
input to every compartment. In addition, it should be noted that the speci-
fication of bounds on the controlled elements implies that a system otherwise
completely controllable must now have a restricted range of controllability in
state space.
The concept of observability of a dynamic system is associated with the
processing of data obtained from measurements on the system. Analogous to the
27
-------
test for controllability, Kalman (1963) has shown that a linear, time-invar-
iant system defined by (5.1) is observable if and only if the n x np matrix
^ = [CT; ATCT; (AT)2 CT; ...; (A1)"-"1^] (5.9)
is of rank n. If there is only one output (p = 1) then a necessary and
sufficient condition for observability is that P^ be nonsingular. The test
for observability may be applied to compartmental models by utilizing tech-
niques that were employed when discussing controllability. Upon inspection of
the definition of observability, it is clear that observability does not
depend upon the controlled input matrix, but rather on the output matrix. A
linear compartment model, in the structural sense, is completely observable if
at least one output is reachable from every compartment. It is also necessary
that no outputs be associated with only the state components which form one or
more sink, all influencing the same outputs (Cobelli and Romanin-Jacur, 1976).
Definition 5.5
A subsystem of compartments which receives input from the remainder of
the system, but has no transfers out of the subsystem from any of its compart-
ments, is called a trap. A trap may also be called a closed subsystem.
It is important to mention that the conditions given for controllability
and observability are not sufficient, that is, the preceding conditions
neither guarantee the identifiability of A, nor are the elements of A_ in the
uncontrollable and unobservable part of tFe system necessarily unidentifiable.
Counter-examples which clarify these points are considered by Distefano
(1977). Bellman and Astrom (1970) emphasize the difficulty of the problem of
identifiability and since that time, only necessary conditions have been re-
lated in the literature. Cobelli and Romanin-Jacur (1976) give criteria for
sufficient conditions, but a close study (Delforge, T977) shows by counter-
example that the proposed criteria are necessary, but not sufficient. The
problem reduces to a recurring mathematical question which asks how many
solutions a set of nonlinear equations has.
Let A be unidentifiable. Then it may be possible to formulate an experi-
ment to show that a subset of compartmental levels can be controllable even
when the system fails to meet controllability and observability requirements
(Johnson, 1976). For example, consider a system in which certain states are,
for biological reasons, not accessable and their evaluation, according to the
structural identifiability tests, is not possible. If a feasible experimental
condition can be obtained by adding a new input, thereby yielding necessary
conditions, the system may be considered identifiable (Milanese, 1976). Also,
it has been suggested that the simplest way to deal with a multicompartment
system is to ignore closed subsystems if possible because any or all of the
traps can be deleted without altering the behavior of the rest of the system
(Thron, 1972). In general, this is possible with regard to measurements on
compartments of completely open subsystems; and most systems utilizing the
direct method for identification of the A matrix would fall into this cate-
gory. However, certain observations, such as the kinetics of urinary excre-
tion, require the use of an indirect measurement and hence the accumulation
pools may not be excluded. Necessary conditions for identifiability are
therefore derived for the indirect algorithm,
28
-------
INDIRECT METHOD
The necessary conditions given for controllability and observability of a
system for the direct identification algorithm imply that any model of a real
system with one or more compartments acting as stores for the system may not
be identifiable. With the indirect scheme, the observer states y_ form a set
of flows disjoint from the internal couplings, each accumulation compartment
being a singular trap. In this section, criteria are investigated for con-
trollability and observability related to compartmental structure of a system
having a accumulation pools (1 < n) from which all observations are to be
made. The local model of (3.1) is again incorporated along with the linear
differential observer for the accumulation pools:
(5.10a)
which can be combined to give
A 0
K 0
from which the total system behavior may be explored,
ment vector z as follows:
z = [0 I]
(5.10b)
Now define the measure-
(5.10c)
Measurements are to be taken from each of the n observer states. The system
described by (5.10) is controllable if and only if the partitioned matrix
B.
0
L- -"
AB
KB
i- '
A2B
KAB
A3B
KA2B
KAn"2B
(5.11)
spans the (a + n) state space on which the system is defined, where the two
n x nrn submatrices P, and j> are called the upper and lower partitions, re-
spectively, of the controllability matrix. Now if P. is of rank 2n, then P,
and F> must both be of full rank. Note that the upper partition of the con-
trollability matrix is the same test matrix as that for the direct method (see
equation (5.8)). Thus, as with the direct method, if there is no path from
any input to compartment i, then the i-th row of the controllability matrix is
identically null, as no input enters compartment i. Next consider the lower
partition of the matrix P_. . If, in addition to the aforementioned condition,
there is no accumulation compartment that can be reached from a given compart-
ment i, then every element of the i-th row of the controllability matrix is
29
-------
null as no pool is influenced by the state variable associated with the i-th
compartment.
For a clearer understanding of the observer portion of the controllabili-
ty matrix, consider the dynamic states y_ to be continuous observations rather
than accumulation pools. Then from Kalman (1963) and the definition of out-
put controllability, the system is said to be completely output controllable
if and only if the composite matrix
P = [KB; KAB; KA2B; ...; KA^B; D] (5.12)
\J\f ~
is of the same rank as the row order of J<, where D_ is the connectivity matrix
which postulates a direct transmission from u^ to y_. The existence of a direct
input into any accumulation pool is not assumed to exist, that is, D^ = 0_. In
fact, if the same compartment is available for insertion and observation of
tracer, identification requires that the associated flow matrix be irreducible
(Smith, 1976), whereas the flow matrix for the indirect identification scheme
is reducible. Thus, tracer must be injected into the internal system ^ and
not into the observer y_. If the rank of the £ x nm matrix
£2 = [£; KB; KAB; KA2]}; ...; KA_n~V| (5.13)
is of rank £, then the rank of the & x nm matrix P given by (5.12) is neces-
sarily £. Notice that state controllability is neither necessary nor suffi-
cient for output controllability.
The discussion of indirect controllability to this point has relied upon
the knowledge of the rank of P. which is not customarily known when discuss-
I \+
ing identifiability. Conditions must be placed upon the submatrices P, and F>
such that the indirect controllability matrix will be of rank £ + n, giving a
necessary condition for identifiability.
Definition 5.6
A system is said to be reduced controllable if and only if the rank of
the matrix P. where
4 = [13; AB; A2B; ...; A1^1]}] (5.14)
is n. The smallest possible integer k for which P^ is of rank n will be
called the reduced order. Sufficient conditions for the controllability
matrix P. to be of the desired rank follow.
ic
Definition 5.7^
Two m x n matrices are defined to be equivalent if and only if one may be
transformed into the other by nonsingluar elementary row and column transfor-
mations. If A_ and B^ are equivalent, we write
A i B (5.15)
30
-------
Theorem 5.1
A system given by (5.1) in which measurements are obtained from £ accumu-
lation pools is controllable if the test matrix P_. (5.11) of the system is
such that 1C
(i) The upper partition P, is a reduced controllability matrix P.
with 2k < n-1, where n is the rank of A and k is the reduced
order of the controllability matrix, and
(ii) CAkP, is of rank £.
Proof: "*
Rewrite the indirect controllability matrix as
B; AB; ...;
Ak+1B; ...;
0; KB; ...; KAk"]B; KAkB; ...; KAn"2B
(5.16)
Now if the matrix P. is reduced controllable, then P, is of rank n where
^ = [B.; AB; ...; AkB] . (5.17)
Since P. is of rank n, then P, is of rank n and
= [I; 0] (5.18)
where I_ is the n-th order identity matrix. Next multiply the matrix F> by
k '
CA , yielding
P£ = [CAkB; CAk+1B; ; CA2kll-
If 2k < n-1, then 2k is at most n-2 so that if P£ is of rank £, then
PJ = [CAkB; CAk+1 B; . . . ; CAn~Vl
(5-19)
(5.20)
is of rank I. (Recall that the rank of a product is less than or equal to the
rank of either matrix. Therefore, the matrix C^ is at most of rank £.) Notice
that F> can be rewritten as
(5.21)
P Pn P B P A B . D"
_O ~ Ly.5 IPs ... 9 ^M J5, - - - 5 £b.
and if P! is of rank n, then
P_2 = [0;
(5.22)
where J^1 is the 5,-th order identity matrix. Therefore
l o
o 1
= I"
(5.23)
31
-------
where _!_" is an (n+£) x (n+£) identity matrix. Thus, the indirect controll-
ability matrix P_-c is of rank (n+£) which is a necessary and sufficient con-
dition for controllability.
Two special cases of the indirect controllability matrix are to be inves-
tigated.
Case 1. If B^ and CB^ are nonsingular, then the matrix P_. is of rank a + n.
This result follows because P. can be regularly partitioned into a triangular
block matrix
p.
E
1C 0 CB
(5.24)
which is nonsingular if and only if each diagonal box is nonsingluar (Parker
and Eaves, I960).
Case 2. If C^ is the n x n accumulative flow matrix of (5.1) which is used in
the indirect algorithm, C_ is nonsingular because all flows are on the diagon-
al. Now if the matrix A_ is multiplied on either side by a nonsingular matrix,
the product has the same rank of A. That is, the rank of CA is n. It
k k
follows that the rank of CA is n and if P. is of rank n, then CA P, is of
k
rank n. In theorem (5.1), the matrix CA FY had to be of the rank C_. There-
fore, in the special case where there is an accumulation pool for each com-
partment, this condition is automatically satisfied.
Theorem 5.2
For indirect identification, there must be at least two inputs to the
compartmental system.
Proof:
The number of accumulation pools within a system is assumed to be less
than or equal to the number of compartments in the system. That is, a <_n.
In addition, the number of rows in P. must be at most equal to the number of
1 C
compartments n times the number of inputs m, or
a + n <_ nm. (5.25)
Now divide both sides of (5.25) by n, yielding
A + 1 < m, (5.26)
n
but since a <_ n, then ^ 5 1 and unless £ = 0 (which is the direct method),
m £ 2 because there can be no partial inputs.
Measurements of the accumulation pools are to be obtained for identifica-
tion of the system of equations (5.10). In order for the system to be
able in the sense of Kalman (1963), it is necessary and sufficient that the
32
-------
test matrix
T T T Tn~2 f
0; K; A'K ; ...; A' K1
I1; 0; 0; ...; 0
is of rank n + £ where _!_' is the £-th order identity matrix.
submatrix P- is of rank n, then
^3 = CO; 1]
where I_ is the n-th order identity matrix and hence,
(5.27)
Note that if the
(5.28)
-3
0 1
r o
1 o
o r
(5.29)
Thus, for P. to be of
giving an equivalent (n+Jl) x (r\+i) identity matrix.
the desired rank, P_~ must be of rank n. By considering the states y_ as con-
tinuous observations of the internal states x_, the rtecessary conditions pre-
sented by Cobelli and Romanin-Jacur (1976) for observability for the direct
algorithm may be applied to the indirect scheme. In the structural sense, the
indirect model is completely observable if at least one flow to an accumula-
tion compartment is reachable from every internal compartment. In addition,
no accumulative flows may be associated with only the internal state variables
which form one or more traps, all influencing the same pool.
COMBINED DIRECT-INDIRECT
Applications of the direct and indirect algorithms may tend to overlap in
certain systems. An introduced (labeled) substance may metabolize and accumu-
late, along with its metabolites, in one or more pools and also be directly
measureable as an output from compartments within the system itself. In other
systems, the outputs available for data collection might be the sum of obser-
vations from a compartment and its accumulation pool.
The method in which the observations are made is not involved in control-
lability which is concerned with the flow of the inputs, so the indirect con-
trollability conditions suffice for the combined measurement technique.
Observability of the two types of combined systems will be discussed in
this section. Requirements for observability are straightforward. First con-
sider the system of (5.10) where the measurement vector is of the form
z =
C_ ^
0 I
(5.30)
Measurements are to be collected from both system compartments and accumula-
tion pools. The system is observable if and only if the test matrix
33
-------
0; C_T
J^
ATCT; ..
0.
0; KT; ATKT; ...; (AT)n'2KT
0
(5.31)
is of the rank p + a where _!_ is an fc-th order identity matrix. Define the
upper partition to be
= TO- P
L-'
(5.32)
and notice that F> is the observability matrix P for the direct algorithm
O 0
while P^ is the observability matrix P_. for the indirect scheme for identifi-
cation. The system is completely observable if the conditions for direct
observability are satisfied because F> would be of rank p, making the system
I1 o
0 I
0
(5.33)
of rank p where V and J^ are p-th order and £-th order identity matrices,
respectively. Also, if there are as many (or more) accumulation pools as
directly measureable states (p), and the conditions for indirect observa-
bility hold, then the system is of the form
1' 0 P|
010
and the desired rank is obtained. In addition, if the two test matrices
(direct and indirect) combine to form a matrix P_7 of rank p, the system
becomes
(5.34)
P_ L7
I 0
E
r o
o i
(5.35)
and the necessary and sufficient conditions for identifiability are fulfilled.
Consider next a system which has as its output the sum of observations.
In this combined case, the measurement vector for the system is
L = EC 1]
(5.36)
and the system is said to be observable in the sense of Kalman (1963) if and
only if the test matrix
Ps = [CT; ATCT; ...; (A1)"-"1!:1] + [0; KT; ATKT; ...; (AT)n-2KT] (5.37)
is of rank p. Thus, the system is observable if the sum of the observability
matrices for the direct and indirect algorithms is of rank p. However, two
34
-------
matrices of different orders cannot be added, so in order for the test matrix
P^ to be defined, the number of direct observations, p, must equal the number
of accumulation pools, £.
35
-------
SECTION 6
CONSTANT INFUSION
A constant infusion input has occasionally been used in tracer methods to
label compartments to a high specific activity. The constant influx of tracer
was probably first incorporated by Hevesey and Hahn (1940) to obtain high con-
32
centrations of P-phosphate in plasma. Hevesey and Hahn, however, were not
interested in analyzing the data to determine model parameters. Constant in-
fusion is also used when there are problems with incomplete and non-instantan-
eous mixing within compartments, usually due to spatial homogeneities. Ship-
ley and Clark (1972) present textbook cases in which the method and its usage
are discussed. Assuming that the chronic infusion of tracer does not disturb
the behavior of the recipient compartments, then the constant input level of
Table 1 gives rise to a steady-state defined by x^.
The constant infusion method as discussed in Section 3 is based upon the
compartment model defined by (3.1) as
x = Ax + Bu_. (6.1)
Linearly independent inputs generate an ensemble of linearly independent
steady-state samples of the system state x^ from which the coefficients of the
matrix of rate coefficients may be determined.
When in steady-state, the system of (6.1) can be described by
0 = Ax_+Bu. (6.2)
Recall that _x is an n-vector of compartmental states, u_ is an m-vector of
specified inputs, and J3 is an n x m matrix where
n
and £ b.. = 1.
i=l 1J
Let u,, u_2, . .., u^ span Rm. Assuming that the compartmental system is open
'1
so that A' exists, then £r Xg» » 4, span at most an m-dimensional subspace
A x = -IB u k=l, ..., m. (6.3)
of Rn where
36
-------
Now set
-By. = w. k=l, ..., m (6.4)
which results in
AX = W (6.5)
where the n x m matrices
X = [xr x2, ..., xj (6.6)
and
W = [wr w2, ..., wj (6.7)
are both of rank m. Taking the transpose of each side of equation (6.5) gives
X1 AT = WT. (6.8)
Next define the n-vector a. to be the k-th column of A_ such that
AT = [a_r a_2, ..., aj (6.9)
The vectors a. define the unknown system parameters, i.e., the rate coeffic-
ients of the compartment model. Thus,
X1 a^ = z^ k=l, ..., n (6.10)
where z. is the m-vector of the k-th column of
Z = WT (6.11)
where
Z = [z_r z2, ..., z^] (6.12)
and rank (I) = m. It should be noted that the vectors a. and z. of (6.10) are
completely determined by the compartmental inputs. That is, a. represents the
k-th row of A_ giving the endogenous inputs to the compartment and z. defines
the possibility of an exogenous input to compartment k.
Identifiability related to the system structure for a multi-input/multi-
output system where all states can be measured directly is to be considered
first. The system is structurally identifiable when the number of parameters
that can be estimated is not smaller than the number of known system parame-
ters in the algebraic equation (6.5).
Only compartment models which are open to their environment are consid-
ered in this chapter because closed systems with constant rate inputs accumu-
late material at a constant rate, never reaching steady-state. This condition
does not prohibit the closure of an individual compartment or groups of com-
37
-------
partments as long as they do not form a decoupled subsystem. The A-matrix for
an open system is characterized by a nonzero determinant. This result implies
that in the special case when the inputs u^ (k=l, ..., n) span Rn, the resul-
tant steady-states x^ (k=l, ..., n) also span Rn and Xf1 exists. Hence, the
matrix equation (6.5) can be solved by multiplying on the right by Xf1, giving
A = WX_ . Input access to every compartment in the system and the ability to
measure all compartmental steady-states then imply structural identifiability.
These sufficient conditions are only rarely met in practice, requiring for
biochemical systems the organic synthesis of tracer material for all chemical
species present in the system. For ecological systems, only a few natural
input pathways to the system exist with many compartments being essentially
inaccessable. Therefore, the remainder of this chapter concerns the develop-
ment of conditions which relate identifiability to system structure with re-
stricted compartmental access.
Theorem 6.1
A compartmental system is identifiable only if there exists at least one
path from.some constant input to every compartment.
Proof:
By contradiction, if there is no path from any input to one or more com-
partments, then for all k, equation (6.3) can be partitioned as
A A
-1 -2
0 A,
o _
Hence, the following relationships arise:
A
But A'1
0
(6.13)
-2
(6.14)
exists, that is, there are no closed subsystems, and the vector x^2 is
identically null. Thus, regardless of the values of the elements comprising
A,, the model parameters contained in A- and A« cannot be identified by sol-
ving equations of the form (6.10). Thus, structural identifiability depends
on the ability of inputs to reach every compartment.
Theorem 6.2
A system is identifiable by the constant infusion method only if at least
one outlet to the environment is reachable from every compartment along any
path.
Proof:
By contradiction, if no system outlets can be reached from compartment m,
then this compartment and all compartments coupled to it continuously^accumu-
late material without reaching a steady-state. This, of course, prohibits
identification using (6.5).
38
-------
The preceding concepts of structural identiflability based upon the data
derived from constant infusion experiments parallel those already known for
the single dose method (see Section 5). However, the solution of the identi-
fication problem of a constant input applied to a system is further con-
strained by the static nature of the resultant steady-state. Indeed, as will
be shown, the dynamic systems concept of observability plays no role in sol-
ving this algebraic problem.
Theorem 6.3
For m exogenous inputs u^ to an identifiable compartmental system, there
can be at most (m-1) endogenous inputs to any single compartment k.
Proof:
For m inputs into a system, equation (6.10) produces m independent equa-
tions in the n unknown parameters of a, . But at most m unknowns can be deter-
mined from the solution of m independent equations. Recall that the k-th
entry in the n-vector a. represents the turnover of the k-th compartment. If
the compartment turnover is to be estimated, then at most m-1 of the endogen-
ous inputs to the k-th compartment (the remaining entries in a. ) can be iden-
tified. Thus, for m inputs to a system, there can be at most m-1 endogenous
inputs to a single compartment.
The (n-m) compartments r without exogenous inputs produce m dependent
equations of the form
XT ar = 0 (6,16)
where a is defined as in equation (6.10). Therefore, (n-m) additional equa-
tions are needed for the solution of (6.5). The necessary equations may be
obtained from each of the (n-m) compartments a in one of the following ways:
1. Turnover of x must be known from an independent measurement, or from
a. pnJLonJi knowldege of the compartmental turnover, giving the r-th
entry in a .
2. Flow out of the compartment to the environment must be determined for
an open compartment.
3. A closed compartment offers an additional independent equation given
by the null column sum
£ air = 0 (6.17)
where a. (1=1, ..., n) define the components of a .
Consider the three compartment system of Figure 9 which gives an A-matrix
of unknown system parameters of the form
39
-------
A =
"an
a21
.0
0
a22
&OO
a!3~
0
av*-
(6.18)
1
x
1
1
V
x2
5
Figure 9. Hypothetical compartment system.
Assume the two exogenous inputs u, and u? drive the system to a steady-state,
resulting in the following equation of the form of (6.12):
21
22
!3
i32j
"V
X2
-XT-
+
1 0
0 1
_0 0_
(6.19)
Performing m=2 experiments with two linearly independent constant inputs gen-
erates the ensemble of steady-states given by (6.5) as
!3
21
22
xn Xi2
X21 X22
'11
J21
'12
22
Lo
Taking the transpose of each side yields
xll X21 X31
X12 X22 X32
21
-a!3
33-
un U2i
U12 U32
(6.20)
(6.21)
40
-------
which is the matrix form of the following three sets of equations in the six
unknown rate coefficients:
A12
for the first column of A,
X31 a!3 "ull
X32 a!3 = "U12
(6.22)
All 21
X12 a21
for the second column, and
X21 a32
X22 a32
X21 a22 "U21
X22 a22 = "U22
X31 a33 = °
X32 333 = °
(6.23)
(6.24)
for the third column.
The system of Figure 9 may be tested for identifiability by inspecting
the necessary conditions provided in this section. There is a path from some
input to every compartment and one outlet to the environment is reachable from
every compartment, thus satisfying Theorems 6.1 and 6.2. The system has two
exogenous inputs and at most one (m-1) endogenous input to a single compart-
ment. Notice that the sets of equations (6.22) and (6.23) have a unique solu-
tion if and only if
ll X31
12 X32
0
and
ll 21
X12 X22
0.
(6.25)
Furthermore, equation (6.24) has a nontrivial solution if and only if
X21 X31
X22 X32
= 0.
These conditions are required for identifiability and are obtained when the X_-
matrix of steady-states has row rank of m=2 and the individual state vectors
can lie anywhere in the m-dimensional subspace of the state space Rn. A suf-
ficient condition for these constraints on the X-matrix of steady-states is
provided by the controllability of (A_, Bj.
Equation (6.22), (6.23) and (6.24) can provide five independent equations
in the six unknown entries in the A-matrix. Thus, with this example, it is
clear that knowledge of the turnover a.,- of compartment 3, or the exogenous
flow (a33)x3 from compartment 3 provides the identification of the systems
model. It is also clear that closure of compartment 3 gives
41
-------
'13
+ a
33
= 0,
(6.27)
which along with equations (6.22) and (6.23) and one equation drawn from
(6.24) defines six independent equations in the six unknown entires in A.
If the compartment(s) without an exogenous input is (are) not closed, it
may be difficult, if not impossible, to identify the compartment model from
steady-state measurements alone. In such cases, it may be useful to define
the maximum identifiable sub-model of the system. The identifiable sub-model
parameter matrix A consists of all transfers associated only with compart-
iTiQ A
ments which are either closed or have exogenous inputs. The dimension of the
sub-model is (m+c) where c is equal to the number of compartments without an
input or a flow to the environment. For example, consider the system of Fig-
ure 9 with an additional flow from compartment 3 to the environment. The
matrix of rate, coefficients for maximum identifiable sub-model is
A,
max
21
22
(6.28)
The identifi ability of this sub-model depends only on the measurement of
steady-states. In order to identify the full system model an experimental
design which incorporates material flow measurements must be implemented.
As a practical matter, it is often true that the closure of a compartment
indicates its inaccessability with respect to exogenous inputs. Thus, the
case when (m+c) = n is important in the applications of the theory of struc-
tural identifi ability. For this case the following theorem gives sufficient
conditions for identifiability.
Theorem 6.4
An open compartment model with m inputs is identifiable provided:
1. The pair (A_, B) is controllable.
. _, .
2. A maximum of Tm-1 ) endogenous inputs enter each of the n compart-
ments.
Proof:
IlltHlUb .
3. The (n-m) compartments without exogenous inputs are closed.
Because (A_, J8) is controllable, the row rank of Y^ is m, and
,,T
k=l , 2,
m
(6.29)
has a unique solution for the m unknowns in the n-vectors a^. Also, the
controllability of (A_, Bj along with the fact that A is open implies that
k=l, 2,
(6.30)
where c=n-m, gives c(m-l) equations in cm unknowns. These results follow from
the fact that controllability implies an input to (6.1) can produce an arbi-
trary steady-state which is not constrained to lie in a proper subspace of R .
Also, since A is open without closed proper subsystems, the resultant coeffic-
42
-------
lent matrix [see (6.24)] of (6.30) is of rank (m-1). So, equations (6.29) and
(6.30) give
m + c(m-l) = mn-c
(6.31)
independent equations in nm unknowns. However, c compartments are closed,
yielding c independent equations for the columns of A. Thus, the solution of
the resultant set of mn linearly independent equations gives the desired nm
rate coefficients comprising the A-matrix.
Now consider the case in which some states cannot be measured directly.
It will be shown by example that it is impossible to identify all unknown
parameters unless there is access to every compartment. Assume that steady-
state measurements can be taken only from compartments 1 and 2 of Figure 9.
The ensemble measurement matrix Y_ may be given by
1 = C.X.
where C is defined by equation (5.2), resulting in
(6.32)
y
11
'21
'12
'22
=
1 0 0
0 1 0
" Xll X12~
X21 X22
. Xo-i X~?j
(6.33)
substituting the measured values Y_ into the system of equation (6.20) gives
the system
'11
21
ln
'13
33
32
'22
'22
-u
-u
0
-u
-u
11
21
12
(6.34)
22
with the additional equation a,~ + a_3 = 0 provided by equation (6.24). No-
tice that the result is a set of six linearly independent equations with eight
unknowns since the values of x--, and
cannot be measured.
Although the parameters associated with compartments which are immeasur-
able cannot be estimated, again, a maximum sub-model can be identified. The
identifiable sub-model flow matrix is composed of those parameters associated
only with those compartments which can be measured directly and satisfy con-
ditions given when all states are measurable. Thus, for Figure 9, the maximum
identifiable sub-model flow matrix A is that of (6.25).
max
Controllability of the system (6.1) allows the exogenous inputs to pro-
duce arbitrary steady-state values anywhere in state space resulting in
43
-------
sufficient, but not necessary, conditions for the structural identifiability
of the system when direct measurements can be obtained from every compartment.
There is, however, no need to investigate the observability of the system as
this concept plays no role in solving the problems of constant infusion iden-
tifiability.
Although the conclusions on identifiability for the constant infusion
method are somewhat limited by system structure and the lack of necessary con-
ditions, it should be pointed out that the prevalent notion that "identifica-
tion of dynamic parameters on the basis of measurements is only possible if
the measurements are taken when the system is in a transient state" (Graupe,
1976, p. 7) is erroneous, because complete sets of parameters can be identi-
fied when a system is in steady-state under conditions which are often met in
practice.
44
-------
SECTION 7
EXPERIMENTAL DESIGN
The experimental design for microcosms receiving toxic material inputs
developed by this research project focuses upon deciding what can be measured
and what models support the measured data. The basic design requires pre-
experiment measurements in order to determine the parameters of the Mulholland-
Gowdy sampling theory for system identification and modeling. Some idea of
the time scale of the system must be obtained for a useful definition of the
flows and material transfers between compartments. Of course, the number of
compartments in the system should be discerned as a first step.
The Mulholland-Gowdy sampling theory requires estimates for the maximum
compartmental turnover and the minimum exogenous transfer rate. The maximum
turnover describes how fast the toxic material moves from one compartment to
another, while the minimum exogenous turnover bounds the toxicant degradation
rate for each compartment. Pre-experiment measurements are also needed to
establish the precision of the measured data, e.g., the number of significant
digits in the data output of the laboratory measuring equipment. The periodic
sampling interval T for time series state data collection resulting from a
basic experiment, as described by Table l(e), is then defined by (3.4) and
(3.6).
Many diagnostic inputs, of the form of basic experiments (Table 1), are
required to perturb the laboratory system in order to elicit the transient
state response required for data collection and system identification. Thus,
the basis for the experimental design is a diagnostic input of the form shown
in Figure 10. It is assumed the laboratory system is of a flow-through type,
open to its environment, so that a constant inflow of toxic material u
applied to the natural microcosm inputs establishes a steady-state. This
initial system steady-state is observed by noting the constant levels of
toxin within each compartment. The inflow u is a chronic level of toxin
carefully chosen to have minimal effects on the biota. The time required to
establish the initial steady-state x response to u is approximately equal to
six (6) time constants, obtained by multiplying the reciprocal of the minimum
exogenous turnover by six (6). Assuming the flow-through process takes place
at a rate greater than degradation, the settling time estimate produces 6/6,
which essentially accounts for the exponential behavior of degradation in the
transient state. At time T-| the input u is reduced by Au, where
u1 = UQ - Au, (7.1)
so that the microcosm seeks a new steady-state x,. In order to identify a
45
-------
'u(t)
T.
'N+l
Figure 10. Diagnostic input for experimental design
xi
x(t)
1 I 4
| . ,
I | | I I
'ill'
' ' ' ' I
Figure 11. Response to diagnostic input for experimental design
46
-------
linear model for the laboratory ecosystem, n samples of the transient state
must be resolved between the steady-states x and x,. Thus, the input change
Au must be large enough to produce this amplitude resolution. For linear com-
partment models, equation (3.2) prescribes that the transient state x(t) is
bounded by
x(t) > x] + (xo-x.,)e-at. (7.2)
Linear system theory also gives
(uo/ul} = (xo/xl}' (7'3)
and since
x1 = XQ - Ax (7.4)
it is easily shown that (7.1), (7.3) and (7.4) combine to give
Ax = (x /u ) Au. (7.5)
Assuming a sampling period of T, it is clear that x(mr), the n-th sample taken
after T,, must be different than x,, i.e.,
x(nr) > x] + 10"d (7.6)
within the finite arithmetic precision of the measurements. Thus, equation
(7.2) gives
Axe-anT > l(Td (7.7)
and a lower bound for Au results by substituting (7.5) into (7.7):
Au > (UO/XQ) e'anT 10~d. (7.8)
This process is continued through time T^ resulting in enough measured data to
construct N (local) linear models for ecosystem structure and function. Since
only positive inflows of toxic material are possible
NAu < UQ (7.9)
must hold true. Hence, an upper bound on Au is obtained
Au < u /N. (7.10)
In response to the staircase diagnostic input of Figure 10, at least n
samples of the microcosm state are measured with a sampling interval of T. As
in Figure 11, the time difference (T.+1 - T.) during which the transient state
response is observed must be greater than six (6) time constants, i.e.,
- T. > 6/6 (7.11)
47
-------
to insure the system is in steady-state before the next input transition takes
place.
The Mulholland-Gowdy sampling theory has been used to design an experi-
ment for studying laboratory ecosystems. Experimental design parameters (a,
6, d, etc.) can all be obtained from pre-experiments involving the microcosm
or its isolated components. The result of an experiment based upon the diag-
nostic input of Figure 10 is a data set which is capable of constructing a
series of (local) linear models for the microcosm. Furthermore, the validity
of these linear models can be tested using the measured data. Two tests exist
for model validation. First, for linear models, the input level ratios must
equal their respective steady-state response level ratios:
for all i = 0, 1 , .. . , N. Of course, practical laboratory experiments will
result in statistical estimates for these steady-states, and standard statis-
tical tests can be used to determine the validity of (7.12). The second test
is qualitative, noting that linear compartment models produce monotone tran-
sition states in response to basic experiments in which transfers between
steady-states are observed. If the samples of the microcosm state during any
- T- ) time period do not decrease monotonically from steady-state x. to
then the linear compartment modeling hypothesis is incorrect. Probably
all microcosms will fail this test for linearity when Au is too large. Thus,
local behavior of the system and the ability to describe this behavior by
linear models are believed to be related.
The validity of the experimental design, as based upon the linear com-
partment modeling objective, depends upon the equality of the steady-state
ratios and the monotonic behavior of the transition between steady-states.
When these validation tests fail, the microcosm modeling objective must be
based upon nonlinear system theory.
If a great deal of intuition exists regarding the microcosm processes,
then it may be possible to characterize these processes mathematically within
a nonlinear model of known structure. In terms of the general concepts of
systems theory, it is extremely difficult to identify nonlinear models using
transient state data. The technique presented in Section 3 of this report
applies only to linear systems and because it relies on a general solution for
the linear model, the scheme cannot be generalized to nonlinear identifica-
tion. However, another approach for obtaining the parameters of a nonlinear
model of known structure is possible using only steady-state measurements.
The proposed method of nonlinear steady-state model identification is
best illustrated by example. Consider the very elementary single compartment
system with input u and state x. From pre-experiment observations, scientific
intuition, or theoretical constructs regarding the microcosm processes, the
structure for the model is assumed to be of the form
r\
x = -ax + bx + u (7,13)
where a and b are unknown parameters which characterize the model structure to
the particular ecological processes. For b = 0 the system is linear and the
parameter a is the compartmental turnover. The identification of the param-
eters a and b proceeds as follows. Let a constant input u1 be applied to the
48
-------
system, which results in a steady- state x, being observed. Then another con-
stant input u2 t u.j is applied, and its corresponding steady-state x2 is also
observed. Substituting these steady-states into (7.13), while noting that a
steady-state is defined by the time derivative of x being identically zero,
two simultaneous algebraic equations result:
2
+ bx, + u,
(7.14)
Because the inputs and states are measured, and hence known, equation (7.14)
represents two linear algebraic equations in the unknowns a and b. Further-
more, a unique solution of (7.14) for a and b exists if and only if x1 r4 x2>
Thus, a model of known structure can be identified from steady-state measure-
ments. However, since an infinite number of nonlinear models could fit the
steady-state data, either confidence in the model structure must be present or
some validation scheme for the model is required.
To demonstrate the multiplicity of nonlinear models which can be identi-
fied from a fixed steady-state data set, consider the model structure:
x = -ex + dx + u (7.15)
where c and d are unknown parameters and k > 3 is any integer- The same
steady-state data set used for identifying (7.13) now gives
0 = -ex, + dx^ + u,
(7.16)
which for each fixed k > 3 yields a unique solution as long as x-j ? x,,. Some
knowledge of the ecological process under study and its theoretical formula-
tion could help to decide which model structure best represents the system.
Also, each model structure could be compared with time series data for the
transient response of the system in order to decide the best model fit to the
data.
The scheme presented for nonlinear identification of the parameters of
models of known structure generalizes. For (7.13) the model consists of a
single equation with two unknown parameters which are obtained from two steady
state measurements. In general, a model consisting of n compartments and M
unknown parameters requires N steady-states measurements, where
nN = M. (7.17)
Thus, in Figure 10, the N levels in the diagnostic input provide a data set
capable of determining M = nN unknown parameters in a nonlinear model of known
structure which represents a microcosm of n compartments. At time T^+1 a
large (positive) step input is again applied to the system to re-establish the
initial state. The time series which results from this step input is of use
49
-------
in validation of the assumed nonlinear model structure.
To complete the nonlinear analysis, it should be noted that the N (local)
linear models which result from the diagnostic staircase input can be used to
compute the material flows between each compartment as they depend upon the
compartment levels. Such tables of global measurements of ecological struc-
ture and function have been used by Smith (1970) to identify the model struc-
ture for an ecosystem. Thus, it may be possible to use the local models to
validate the nonlinear model structure, and by doing so provide an important
connection between the linear and nonlinear model identification methods.
A pre-experiment to determine the parameters of the Mulholland-Gowdy sam-
pling theory has been run at the Environmental Research Laboratory in Athens,
Georgia, for a microcosm consisting of water, algae, snails and sediment com-
partments. A slug of radiolabelled methoxychlor was introduced into the water
of the laboratory ecosystem. Samples of the toxin and radioactivity were
obtained in order to determine estimates for degradation and compartmental
transfers. These estimates define a and 6, the maximum turnover and minimum
exogenous turnover. By analyzing the methoxychlor samples, some idea of the
arithmetic precision of the data is obtained. All of these estimates for the
sampling theory parameters determine an experimental design, based upon the
diagnostic input of Figure 10, for the particular microcosm under study. It
is hoped the experimental design will be tested in the laboratory sometime in
the future.
50
-------
SECTION 8
COMPUTER PROGRAM
A digital computer user package implementing the compartment model iden-
tification algorithm has been developed. Programmed in standard FORTRAN, the
code which implements the identification algorithm is available in two forms,
batch and interactive, both of which run efficiently on a standard IBM 370
operating system. The code is independent of any special library routines,
and as such is completely transportable to any IBM System/360 or System/370
series computer. The interactive program can be implemented on any IBM
System/370 time share operating system (TSO).
The software consists of a main program and one subroutine, EIVECT, the
flow diagrams of which are shown in Figures 12 and 13.
The program begins by reading the number N of compartments in the bio-
logical system being analyzed, the number K of measurements obtained from a
single experiment, and the pre-calculated sample period T. Two flags are then
set which indicate row or column data input and printing options, respective-
ly. The program then moves immediately to a loop that reads and stores the
sample data array in matrix form, i.e., Xjt). Since N+l discrete measurements
of the state are necessary and sufficient for the determination of the state
transition matrix $_, the program terminates for K < N+l. If K = N+l, two
matrices of coefficients, X, and )L, are formed from the N+l state samples as
described by (3.12) and (3,13). The solution for the state transition matrix
<£ is acquired by inverting matrix X, of (3.12). Gaussian elimination with
partial pivoting on a matrix scale is used. Upon computing the N x N matrix
subroutine EIVECT is called for the identification of the matrix of rate
coefficients A_. The indices I and L of the main flow diagram are utilized in
the formation of a loop in which several A matrices may be computed for a
comparative analysis if K > N+l.
Subroutine EIVECT solves (3.14) for the matrix of transfer constraints A.
The EISPACK (IBM/360 version by Smith, et al., 1976) is utilized to obtain the
eigenvalues w-j (i = 1, ..., N) and the eigenvectors of <^. The EISPACK uses
QR transformation and decomposition, an explanation of which is given by Wil-
kinson (1965) and Wilkinson and Reinsch (1971).
The matrix $ can be reduced to a diagonal matrix g in which the eigen-
values u>- (i = 1, ..., N) of
-------
( EIVECT J
READ'
PRINT:
CALL
EIGENP
PRINT:
PRINT:
AJJ 0
Aii Xi
A TAT'1
PRINT
A
( RETURN }
\.
Figure 12. Flow chart for computer algorithm.
52
-------
READ:
N, K,X(t)
II 1 + 1
X2 X
CALL
EIVECT
I 1*1
Figure 13. Flow chart for subroutine EIVECT.
53
-------
ST AT = [S A]T (8.2)
where S^ is a similarity transformation having the eigenvalues of $_ as its T
column vectors. Equation (8.2) is solved successively for the columns of A ,
again using Gaussian elimination. Note that the column vectors must be lin-
early independent for the existence of S_ . This follows from the assumption
that the eigenvalues of £ are distinct.
It should be mentioned that EIVECT derives the complex matrix A_ for a
real matrix <£ using the complex principal value natural logarithm in (8.1).
Also, all manipulations are in double precision. The input format for N, K
and the two flags is 14, T is F10.0, and the state data is 8F10.0. The output
is available in full double precision, but the current version of the program
truncates to seven digits as given by 8D14.7 format.
The code, as submitted to the project officer, is well documented by
comment cards, and the particular deck submitted will include data input cards
for Smith's (1970) ecosystem model. An output from this deck, as run on the
Oklahoma State University computer is included in the user package. A com-
plete listing of the code which implements the identification algorithm is
also in the user package.
The user package, as described, fulfills the grant obligations. However,
because the interactive version of the program was of use to the researchers
during the grant period, a listing of this code fs also included. This list-
ing indicates the changes which have been made to the batch program in order
to implement it as an interactive program. Depending on the computer opera-
ting system, it is possible to load the batch program version as a data set
and then edit it to conform to the interactive mode. This was the procedure
used successfully by the researchers at Oklahoma State University.
54
-------
REFERENCES
Bellman, R. and K. J. Astrom. 1970. On Structural Identiflability. Math.
Biosci. 7:329.
Bowden, R. 1973. The Theory of Parametric Identification. Econometrics.
41:152.
Brockway, D. L., J. Hill IV, J. R. Maudsley and R. R. Lassiter. 1978. Devel-
opment, Replicability, and Modeling of Naturally Derived Microcosms.
J. Environ. Studies, to appear.
Carlson, A. B. 1975. Communication Systems. Second edition. New York:
McGraw-Hill.
Chen, C. 1970. Introduction to Linear System Theory. New York: Holt, Rine-
hart and Winston.
Cobelli, C. and 6. Romanin-Jacur. 1976. On the Structural Identiflability of
Biological Compartment Systems in a General Input-Output Configuration.
Math. Biosci. 30:139.
Del forge, J. 1977. The Problem of Structural Identifiability of a Linear
Compartmental System: Solved or Not? Math. Biosci. 36:119.
Distefano, J. J. 1977. On the Relationships Between Structural Identifiabil-
ity and the Controllability and Observability Properties. IEEE Trans.
Autom. Control. AC-22:652.
Draggan, S. 1976. Role of Microcosms in Ecological Research. Bioscience.
26:402.
Emanuel, W. R., D.C. West and H. H. Shugart. 1978. Spectral Analysis of
Forest Model Time Series. Ecol. Modelling. 4:313.
Glover, K. and J. C. Willems. 1974. Parametrizations of Linear Dynamical
Systems: Canonical Forms and Identifiability. IEEE Trans. Autom.
Control. AC-19:640.
Graupe, D. 1975. Identification of Systems. New York: Robert E, Krieger.
Hearon, J. I. 1963. Theorems on Linear Systems. Ann. N.Y- Acad. Sci.
108:36.
55
-------
Hevessey, G. and L. Hahn. 1940. Turnover of Lecithin, Cephalin and Sphingo-
myelin. Biol. Meddr. 15:5.
Hildebrand, F. B. 1956. Introduction to Numerical Analysis. New York:
McGraw-Hill.
Johnson, L. E. 1976. Control and Controllability in Compartmental Systems.
Math. Biosci. 30:181.
Kalman, R. E. 1963. Mathematical Description of Linear Dynamical Systems.
SIAM J. Control. 1:153.
Lee, R. C. K. 1964. Optimal Estimation, Identification, and Control. Cam-
bridge, MA: MIT Press.
Lin, P. L. and Y. P. Yu. 1977. Identification of Linear Time-Invariant
Systems Using Exponential Curve Fitting. Proc. IEEE. 65:1508.
Luenberger, D. G. 1964. Observing the State of a Linear System. IEEE Trans.
Mil. Elec. MIL-8:74.
Mulholland, R. J. and C. M. Gowdy. 1977. Theory and Application of the Mea-
surement of Structure and Determination of Function for Laboratory Eco-
systems. J. Theor, Biol. 69:321.
Odum, E. P. 1971. Fundamentals of Ecology. Third edition. Philadelphia:
Saunders.
Ogata, K. 1967. State Space Analysis of Control Systems. Engelwood Cliffs,
NJ: Prentice-Hall.
Platt, T. and K. L. Denman. 1975. Spectral Analysis in Ecology. Ann. Rev.
Ecol. System. 6:189.
Reid, J. G. 1977. Structural Identifiability in Linear Time-Invariant Sys-
tems. IEEE Trans. Autom, Control. AC-22:242.
Shipley, R. A. and R. E. Clark. 1972. Tracer Method for In Vivo Kinetics.
New York: Academic Press.
Shugart, H. H., ed. 1978. Time Series and Ecological Processes. Philadel-
phia: SIAM Publ.
Smith, F. E. 1970. Analysis of Ecosystems. In; Analysis of Temperate For-
est Ecosystems, D. E. Reichle, ed. New York: Springer-Verlag, pp. 7-18.
Smith, W. D. and R. R. Mohler. 1976. Necessary and Sufficient Condition in
the Tracer Determination of Compartmental System Order. J. Theor. Biol.
57:1.
Thron, C. D. 1972. Structure and Kinetic Behavior of Linear Multicompartment
Systems. Bull. Math. Biophys. 34:277.
56
-------
Tse, E. 1973. Information Matrix and Local Identifiability of Parameters.
In: Proc. Joint Autom. Control Conf., Columbus, Ohio. pp. 611-619.
Tse, E. 1978. A Quantitative Measure of Identifiability. IEEE Trans. Sys.
Man Cyber. SMC-8:1.
Van Voris, P., R. V. O'Neill, H. H. Shugart, and W. R. EmanuelI. 1978.
Functional Complexity and Ecosystem Stability: An Experimental Approach.
ORNL/TM-6199, Oak Ridge National Laboratory, Oak Ridge, TN.
Waide, J. B. 1974. A Linear Systems Analysis of the Calcium Cycle in a
Forested Watershed Ecosystem. In: Progress in Theoretical Biology,
R. Rosen and F. M. Shell, ed. New York: Academic Press, pp. 261-345.
Wilkinson, J. H, 1963, The Algebraic Eigenvalue Problem. London: Oxford
University Press.
Wilkinson, J. H. and C. Reinsch. 1971. Linear Algebra. New York: Springer-
Verlag.
57
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/V78-10*>
3. RECIPIENT'S ACCESSION NO.
4. TITLE ANDSUBTITLE
New Sampling Theory for Measuring Ecosystem Structure
5. REPORT DATE
December 1978 issuing date
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
Robert J. Mulhoi land
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
School of Electrical Engineering
Oklahoma State University
Stillwater, OK 74074
10. PROGRAM ELEMENT NO.
1LA760
11. CONTRACT/GRANT NO.
R-805564010
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Research Laboratory - Athens, GA
Office of Research and Development
U.S. Environmental Protection Agency
Athens, GA 30605
13. TYPE OF REPORT AND PERIOD COVERED
Final. 10/77-10/78
14. SPONSORING AGENCY CODE
EPA/600/01
15. SUPPLEMENTARY NOTES
16. ABSTRACT
This research considered the application of systems analysis to the study
of laboratory ecosystems. The work concerned the development of a methodology
which was shown to be useful in the design of laboratory experiments, the pro-
cessing and interpretation of the results of these experiments, the development
of model structures for microcosms based upon experimental data, the testing
of the predictive capabilities of evaluative models, and the analysis of measure-
ment errors in the presence of noise. Development of the experimental design
methodology was based upon a non-Nyquist sampling theory which gave a priori
computable bounds on the sample period for data collection.
Results of the research project included a definition for the bandwidth
of an ecosystem in the frequency domain, an indication of the effect of measure-
ment noise on the sample period for data acquisition, a scheme for reducing
this noise, and a computer code for model identification based upon the experi-
mental design developed. The project also produced the conditions under which
two new schemes for model identification were applicable to the construction
of linear compartment models. These new identification schemes were based
upon system input perturbations and state observations which were often readily
implemented in the laboratory.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Systems analysis
Mathematical models
Environments
62B
68D
72E
3. DISTRIBUTION STATEMENT
Release to Public
19. SECURITY CLASS (ThisReport)
SECURITY CLASS f
Unclassified
21. NO. OF PAGES
66
20. SECURITY CLASS (Thispage)
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
58
'_ U.S. GOVERNMENT PRIMING OFFICE. 1979 -657-060/1556
------- |