-502-
Finding Qualitative Behavior Modes: The Use of Interval Analysis
to Perform Sensitivity, Stability, and Error Analysis
on Dynamic Models
Weldon A. Lodwick
University of Colorado-Denver
Ralph Levine
Michigan State University
ABSTRACT
A method is described and illustrated for explicit incorporation of and
computation with ranges of initial conditions, functions, and parameter values
in dynamic models using interval analysis. This approach is neither a statistical
nor fuzzy set analysis but instead utilizes interval arithmetic which is
particularly well suited for computerization. When a dynamic model is couched
in interval analytic terms, ranges of all possible solutions are generated
allowing not only an analysis of ranges of behavior modes but for sensitivity
and stability analysis to be performed as a natural part of the model. Moreover,
uncertainties such as specification, numerical method (e.g., numerical
integration), and roundoff errors can also be analyzed in conjunction with or
separate from the interval dynamic model.
INTRODUCTION
In analyzing system dynamic models, time is spent in searching for the conditions,
if any, which generate different qualitative behavior modes. There are a number
of reasons for being interested in the possibility of the model generating two
or more behavior modes. In building confidence in models--see, e.g., Forrester
and Senge (1980)--one would like to reject the models if changes in parameters
yield anomalous behavior that appears to be unjustifiable. Secondly, the real
problem that is being modeled may be multifacited so if the system itself is
known to display more than one behavior mode, the system dynamicist must show
that the model can capture these known dynamics. In cases where the system is
very complex and a priori knowledge about the system relatively unknown, a good
model might generate surprise behavior patterns, which actually may be contained
in the data, but was not discovered until the model pointed out its existence
and perhaps its importance.
Another important reason for seeking methods for aiding in discovering different
behavior modes is to bolster the modeler's argument against stressing the
quantitative aspects of the problem, instead of paying more attention to the
relationship between dynamic structure and the qualitative behavior of the system.
In particular, if a modeler can show that under all reasonable ranges of parameter
values, the model behaves in the same manner, and thus is insensitive to parameter
changes, then the client, editor, and/or grant reviewer may feel less inclined
to insist on elaborate aggregated statistical analysis of the past history of
the system as the sole criterion of the model's validity.
-503-
Moreover, in recent years, there is a growing trend in the system dynamic
literature to explore ways of performing sensitivity analysis on continuous
nonlinear dynamic models (Tank-Nielsen 1980; Forrester 1983; Ford, Amlin, and
Backus 1983; Graham and Pugh 1983). Obviously, methods for searching for
different qualitative behavior modes would fall under the general heading of
sensitivity analysis. Forrester (1983) and Graham and Pugh (1983) take an
eigenvalue approach of linearizing the model around important operating points
to assess possible behavior modes and to attempt to understand the role of each
major loop in the model. On the other hand, Ford, Amlin, and Backus (1983) take
amore global statistical approach to parameter sensitivity, which is somewhat
similar to our method which uses interval analysis. In any event, all approaches
to sensitivity analysis have to contend with several problems:
(1) In large complex models, the potential number of combinations of initial
states, parameter values, and shapes of table functions is enormous. In
fact as will be seen, if in a linear model there are n (constant) parameters
whose values are "fuzzy" but contained in a bounded interval, there are 2
combinations. For continuous functions or nonlinear models whose "fuzzy"
domains and ranges are contained in bounded sets, there are, of course,
uncountably infinite combinations.
(2) Sensitivity analysis is frequently performed in situations where one must
contend with noisy input functions, roundoff and numerical (e.g.,
integration) error, as well as specification error.
(3) Attempting the classical strategy of keeping everything constant except one
parameter or initial value may be somewhat misleading due to the interaction
among parameters. Some qualitative behaviors may only show up under certain
combination of parameter conditions.
We are just learning to apply interval analysis to overcome these problems. The
key idea is to locate those aspects of the model which are the most uncertain
or fuzzy and to represent those uncertainties in terms of bands or intervals of
values. For example, in specifying table functions, usually there are a few
points where, from the logic of the situation, there is no uncertainty at all.
The origin or the point (1,1) frequently are set up in multipliers to cut off a
rate equation when inventories are zero or to insure that the multiplier has no
influence under normal baseline conditions. On the other hand, parts of table
functions may be extremely fuzzy. Interval analysis allows us to put bands
around the more uncertain parts of table functions before the system is simulated,
and therefore brings in a form of sensitivity analysis at an early stage of the
model building process.
One of the potential strengths of interval analysis is its ability to handle any
combination of bands of parameter values and to generate an envelope of output
for each requested variable. The information per run about potential behavior
modes is maximal, and, as we will demonstrate in a simple example, it would take
Many more runs with the classical approach to get similar information. In
addition, we will also show how interval analysis can work with combination of
parameter changes, and how those results could differ from changing one parameter
at a time, even when using the extreme values of the parameters in each run.
-504-
INTERVAL ANALYSIS--INTRODUCTION
Interval arithmetic is an algebraic and topological structure (Moore 1979 or
Alefeld and Hertzberger 1983) which uses as its basic element of analysis closed
and bounded real valued intervals A= Ca) ,a5] = {x: ay SX < ay where ay .8),x€real
numbers}. The interval A is regarded as a number in much the same way as
x = atbi is regarded as a number. Interval numbers however not only possess
an algebraic structure (different from real and complex numbers since the
distributive law does not hold but a sub-distributive one) but, since an interval
number is also a Set, possesses a topological structure as well. We will be
most interested in exploiting the algebraic structure as a tool in systems dynamic
(SD). Since real numbers are a subset of interval numbers; that is, A= a =[a,a]
is areal number (an interval of zero width); the methods developed herein are
an extension of the "usual" SD approaches.
The following notation will be used. When the context is clear, the capitalized
letters A, B, C, X, Y, and Z denote interval variables. Otherwise, a superscript
"I" is used; e.g., al = [aj,a9]. The capitalized letters F, G, and H will
denote interval valued function counterparts to the real valued functions f, g,
and h, respectively, and defined in what follows. When the context is not clear,
an interval valued function is denoted as flock),
The four categories of errors to be analyzed are: (i) model-specification (or
simply specification), (ii) data or measurement, (iii) numerical method or
discretization, and (iv) truncation. Truncation error is of two types:
(iva) numerical (finite numerical representation of an infinite numerical
process), and (ivb) roundoff or finite state machine errors. The methods of
interval analysis, which include interval arithmetic, will be used to explicitly
incorporate errors (ii), (iii), and (iv) in dynamic models. However, when
quantitative (or qualitative) bounds are known for (i), it too can be analyzed.
It should be clear that quantitative (or qualitative) error bounds yield
quantitative (or qualitative) solutions where by quantitative bounds is meant
that bounds are known precisely.
INTERVAL ANALYSIS--A SHORT REVIEW
Though interval analysis has its roots in Archimedes' method to derive an
approximation for 1, its formal study in numerical analysis began in 1962 when
R. E. Moore published his Ph.D. thesis. An extensive literature now exists with
over 757 citations in Moore (1979) alone. The interested reader is directed to
the more recent Alefeld and Herzberger (1983) and Moore (1979) expositions. An
extremely abbreviated synopsis is presented here.
Interval arithmetic is defined as follows. Let A= fa, ,a] and B= [b;,b5].
AB = [aytb,,a,tbo] (1)
A-B = [ay-by,ap-bz] (2)
AB = [min{c},max{c}] , (3)
where c = {a,b ,ab5,apb),apbo}. If O¢ 8B, then
-505-
[min{c},max{c}] , (4)
{a1/b,a4/bp,a5/by,a9/bo}.
A/B
where c
Assuming that the algebraic operations can be carried out exactly, each operation
(1)-(4) above has the property that for any xeA, yeB, xey © AoB, where o
denotes one of the four algebraic operations. On the other hand, for any ze AoB,
there exists x € A and y €B such that z= xey. However, this exact
interval arithmetic (EIA) will not satisfy the above remarks when implemented on
a computer. The presence of roundoff error in computer floating point
representation of real numbers and in computer arithmetic cannot be guaranteed
to produce the exact endpoints. By analyzing the accuracy of the floating point
representation and arithmetic operations on a given machine, it is possible to
"round" the machine-computed endpoints and floating point numbers to insure
containment. We no longer have equality as in (1)-(4). The implementation of
interval arithmetic in the presence of roundoff errors of particular machines
is called rounded interval arithmetic (RIA) and FORTRAN-based RIA exist for CDC,
DEC, Honeywell, IBM, and UNIVAC computers (Moore 1979, pp. 14-17) as well as
general ALGOL 60 RIA's (Alefeld and Herzberger 1983, pp. 288-295).
The essential condition of RIA is that AoB < ASB where o is (1)-(4) and 3
indicates the use of RIA. It follows from monotonicity that, if F is a rational
expression in the interval variables Xp oXoseeeeXns then F(X] sXpaeee9Xq) c
F(Xq sXpy 009 Xp) where F indicates the use of RIA in the evaluation of F. Suppose
that f is a real valued function defined on the interval A. Consider next the
problem of finding the exact values f(x), x ¢A on acomputer. First, x might
not be representable exactly and would itself have to be replaced by an interval
containing x. Second, f might (usually does) involve arithmetic operations
leading to roundoff errors. Third, f might contain irrational expressions and
even be transcendental leading to approximation procedures when computed.
The concept of interval extension (or range) of functions is introduced to deal
with these problems. The interval valued function F of n variables XqoXoseoe eX,
is an interval extension of the function f if F(Xq sXpaeeeoXq) = F(Xy aXoaeee ay)
To simplify the exposition, we treat the one variable case, where the extension
to the n variable case is easily accomplished. The range of a real valued
function f(x), x €X, is denoted RF = {f(x): x eX}. It is proved in Moore
(1979, p. 22) that if F is a rational interval function and an interval extension
of f, then Rf CF. Therefore, at least formally, there is a means of dealing
with the “usual” SD models possessing rational functions and table functions.
In general under. some conditions on f, f can be approximated by a rational
function where the bounds of approximation are well known. These bounds are
added (and subtracted) to form an interval extension formulation.
However, serious problems exist in obtaining a "suitable" interval extension
function F for f in the sense that F(X) is as close as possible to RF; i.e.,
such that pLRf-F(X)] < © where p is the measure (absolute value) on the real
numbers, €arbitrarily small. The following example from Alefeld and Herzberger
(1983, p. 25) illustrates the problem.
-506-
Let f(x) = x(1-x) where we have F. 10%) = X- x2, Fo(x) = X(1-X), F(X) =
4 - (X-4)(X-4), and F(X) = 4d - a i? as interval extensions of f(x),
x eX = [0,1]. Rf = [0,4], and applying (1)-(3), we obtain F(X) = [-1,1],
F 9X) = [0,1], F 36%) = [0,4], and Fy (X) = [0,4]. As it turns out, to obtain
a ? onininal® interval extension, centered forms such as Fy (X) above are usually
used (Ratschek and Rokne 1984, pp. 30-92, 103-110). It is assumed in what follows
that by interval extension F, the minimal interval extension is meant.
SENSITIVITY AND STABILITY ANALYSIS
Sensitivity Analysis: Once the dynamic model is couched in interval analytic
terms and arithmetically performed using EIA or RIA, sensitivity analysis can
be done very simply. This is because variations in the parameters, initial
conditions, and even the error bounds themselves can be viewed as intervals. If
all interval variations are zero width intervals, "traditional" sensitivity
analysis is obtained. However, when one or more are intervals of width greater
than zero, the full range of all possible solutions results. Therefore, interval
dynamic models are an extension of "traditional" sensitivity analysis. Actually,
in view of the incorporation of truncation and numerical errors, the computed
solution is usually (always) a larger range of values containing the actual range
of values of the solution.
Stability Analysis: There are two ways which this analysis can proceed,
empirically or directly. Once an equilibrium point has been found, the empirical
method takes variations in the parameters using various interval values of the
parameters and executes the interval version of the dynamic model to see how the
final solution is changed. Thus, the empirical method is a sensitivity analysis
about the equilibrium point and can be carried out using the techniques already
described.
The direct method requires an interval analytic eigenvalue solver. For example,
the power method could be used in the following way. Implement the power method
in the usual fashion with real numbers and usual arithmetic. Once the algorithm
reaches a predetermined tolerance, perform the power method one more time using
RIA and theoretical error bounds where the first inputted value in the interval
power method is an interval of zero width equal to the last iterate of the real
number power method. Then perform one step of the interval power method. If the
QR algorithm is used, an interval version would have to be used. Lastly, if the
dimension of the mabrix Mt is not prohibitively large, then the zeros of
pal )= det (al -vl I 1 (5
can be solved by techniques found in Alefeld and Herzberger (1983, pp. 101-112).
TWO EXAMPLE IMPLEMENTATIONS
Salient features of the interval analytic approach to SD are demonstrated by the
following two simple examples. EIA is used here merely to simplify the analysis.
A Simple, Linear One-Level Discrete Dynamic Model: To illustrate the application
of interval analysis, consider
=~507—
x(k+1) = aex(k) +b, (6)
x(0) = 1.0, a= -1.0, b=4.0; k =0,1,...,5,
whose properties are known (Goldberg, 1958).
This is an extremely simple difference equation which is not often seen or used
in system dynamics. First it should be noted that the model with those specific
values of a and b will generate oscillatory behavior, even though it is of first
order. This points out the qualitative differences between the dynamics of
differential and difference equations. Oscillatory behavior does not emerge in
a first-order differential equation system, but can arise in a system described
by a difference equation, such as the one above. A first-order difference
equation might be useful in describing how people manage systems in which the
only information comes at a sampled interval. For example, the second author
has worked with doctors who get information about their diabetic patients on a
daily basis. A poor doctor might only pay attention to one state variable, such
as amount of sugar in the blood. If the sampling rate is too slow, the doctor
can cause tremendous oscillations in the patient by prescribing more insulin to
counter the rising blood sugar level.
The output of this oscillatory model would be 3.0, 1.0, 3.0, 1.0, etc. How would
one apply interval analysis to this simple model? The first step is to define
the intervals around each of the parameters and the initial value of x.
Suppose that ae A= (a,,a5], be B= Cb, ,bo], and x(0)e X(0) = [x1 (0),x5(0)].
Then, X(k+1) = AxX(k) +B, where X(k) = [xq (k) X5(k)I and x and + are performed
using (3) and (1). Given variations so that ay <0 and by > 0, the interval
model becomes
xy (k+l) = ay #Xo(k) +b, if Xy(k) > 0
ayeXo(k) + by if Xp(k) < 0
ageXy (k) + bo if xz (k) 20
ayaxy(k) + by iF x(k) <0
Let A= [-1.1,-0.9], B = [3.9,4.1], and X(0) = [0.9,1.1]; then (7) yields
Table 1, where the lower bound x(k), upper bound xo(k)s center = 2(x5(k) -
x1 (k)), and width = (X5(k) - x1 (k)) of the generated interval [xq (k),X5(k)]
are listed. The interval Ex, (k) x5(k)] in Table 1 is the envelope of all possible
trajectories generated by the simultaneous variations in the parameters a, b,
and initial value x(0) given by A, B, and X(0). To attempt to duplicate the
information gained in the one implementation of the interval model (7), the
classical approach might be performed six times, taking the extreme value at
each end of the three variations A, B, and X(0) and keeping all other values
constant. This is done and listed in Table 2.
(7)
Xo(k+1) =
-508-
k 1 2 3 4 5 6 wince 24 25
x1 (k) 2.69 0.281 2.0531 -0.3318 1.42257 -1.01149 ... -18.8693 -9.76447
Xo(k) 3.29 1.679 3.8471 2.25221 4.46499 2.81969 ... 12.4222 24.852
Center 2.99 0.98 2.9501 0.96020 2.94378 0.9410 ... -3.22352 7.54587
Width 0.6 1.398 1.794 2.58402 3.04242 3.83118 ... 31.2915 34.6207
TABLE 1--Envelope of Trajectories
The six classical sensitivity analyses for (6) would be to consider the following
substitutions one at a time holding other values at their original values:
a=-l.l, a=-0.9, b=3.9, b=4.1, x(0) =0.9, and x(0) =1.1, which would
be A=[-1.1,-1.1], B= [4.0,4.0], X(0) =£1.0,1.0]; A =[-0.9,-0.9], and so
on when (7) is used. The results are (min and max values underlined):
k 1 2 3 4 § 6 wee 24 25
a= -l.d;
x(k) 2.9 0.81 3.109 0.5801 3.36189 0.30192 ... -7.0069 11.7076
a= -0.9;
x(k) 3.1 1,21 2.911 1.3801 2.75791 1.51788 ... 2.0171 2.18461
b = 3.9;
x(k) 2.9 1.0 2.9 1.0 2.9 1.0 wa 140 2.9
b = 4.1;
x(k) 3.1 1.0 3.1 1.0 asl 1.0 wee 1.0 3.1
x(0) = 0.9;
x(k) 3.9 0.9 3.1 0.9 3.1 0.9 wee 0.9 3.1
x(0) = 1.1;
x(k) 2.9 1.1 2.9 1.1 2.9 11 ve 1,1 2.9
Center 3.0 1.01 3.0045 0.9801 3.0599 0.9099 ... -2.4949 6.9461
Width 0.2 0.4 0.209 0.8 0.604 1.216 ... 9.0240 9.523
TABLE 2--Six Classical Sensitivity Analyses
If simultaneous variations of (6) were executed with the left and right bounds
of A, B, and X(0); i.e., using a=-1.1, b= 3.9, and x= 0.9, and then using
a= -0.9, b= 4.1, and x(0) = 1.1; we obtain the following:
k 1 2 3 4 5 6 wars 24 25
x(k) 2.91 0.699 3.1311 0.4558 3.3986 0.1615 ... -7.5705 12.2275
x(k) 3.11 1.301 2.9291 1.4638 2.7826 1.5957 ... 2.0735 2.2338
Center 3.01 1.0 3.0301 0.9598 3.0908 0.8786 ... -2.7485 7.2307
Width 0.2 0.602 0.202 1.0080 0.6163 1.4342 ... 9.644 9.9937
TABLE 3--Simultaneous Variations
Note that model (6) implementations were unable to capture the full effect of
the interaction of the errors as can be seen by comparing Table 1 to Tables 2 and
-509-
3. To capture the complete interaction of three pairs of variations [ay,a,],
[b,,b)], and [x1 (0),x9(0)] of the constant parameters a, b, and x(0), eight runs
of (6) must be made and in general a linear model with n interval variations of
(constant) parameters requires at least 2" runs to guarantee the capture of the
full effect of interactions of the n variations. [€ requires ‘one run of the
interval model to obtain the same information. When the model has a continuous
function which changes monotonicity over the variations and subsequent generated
values of the state variable, more than 2” runs (in fact an uncountably infinite
number in many cases) would be required to guarantee one had obtained all possible
interactions of the variations, but requires one run of the interval model.
A Predator-Prey Model--The Kaibab Plateau Model (Goodman 1980, pp. 377-388):
The second example presented is closer to home. This model has been used
frequently in the system dynamics literature for teaching purposes and a DYNAMO
version can be found in Goodman (1980). There are three levels in this model:
deer, predators, and food. The inputs are total area (AREA) = 800,000 acres,
average food per deer (AFPD) = 1.0 units, and the removal rate of predator (RF).
The time period starts in 1900 with the deer population (DP) = 4000, predator
population (PP) = 8000, food capacity (FCAP) = 350,000, and available food for
deer (F) = 350,000 units. In 1905 a bounty on predators is given so that RF = 0.2
from 1905 on. The time increment (DT) = 0.1 is used.
The equations of the model are:
Food per deer FPD = F/DP (8)
Food ratio FR = FPD/AFPD (9)
(1) DEER SECTOR
Deer growth rate factor DGRF = DGRFT(FR) Table Function (10)
Deer net growth rate DNGR = DP«DGRF (11)
Deer density DD = DP/AREA (12)
Deer kill ratio DKR = DKRT(DD) Table Function (13)
Deer predator rate DPR = PPxDKR (14)
(2) PREDATOR SECTOR
Predator growth rate factor PGRF = PGRFT(DKR) Table Function (15)
Predator net growth rate PNGR = PPxPGRF (16)
Predator bounty removal PBR = PPxRF (17)
where RF = 0.0 for 1900 until 1905 and RF = 0.2 from 1905 on.
(3) FOOD SECTOR
Food capacity fraction FTFCAP = F/FCAP (18)
Food regeneration rate (years) FRT = FRTT(FTFCAP) Table Function (19)
Growth rate (units) GR = (FCAP-F)/FRT (20)
Food consumption rate per deer FCPD = FCPDT(FR) Table Function (21)
Food consumption (units) FC = DPxFCPD (22)
(4) INCREMENT
Deer population DP: = DP + DT(DNGR-DPR) (23)
Predator population PP: = PP + DT(PNGR-PBR) (24)
Food available F: = F + DT(GR-FC) (25)
The interval analytic model at first glance would appear to be a straightforward
substitution of each variable with interval variables, each algebraic operation
with RIA (or as in our example EJA), and each table function with an interval
table function. However, if this were done, intervals whose widths are too large
-510-
ensue. This is most clearly seen in (17) and (24). Let the interval containing
PBR at period k be denoted [PBR1(k),PBR2(k)] and likewise for the other variables.
Since P >O and RF > 0, using interval multiply (3), (17) becomes:
PBRI(k) = PP1(k)«RF1(k) (26)
PBR2(k) = PP2(k)aRF2(k) . (27)
Using interval add (1) and subtract (2), (24) in a straightforward substitution
becomes:
PP1(k+1) = PP1(k) + DT(PNGRI(k) - PP2(k)«RF2(k)) (28)
PP2(k+1) = PP2(k) + DT(PNGR2(k) - PP1(k)*RF1(k)) . (29)
However, the interval arithmetic is removing a predator bounty kill PP2(k)«RF2(k)
associated with the largest value of the predator population PP2(k) from the
smallest value of the predator population PP1(k). Clearly, PP1(k) does not
experience that removal PP2(k)«RF2(k). It does however experience a maximal
removal of PP1(k)«RF2(k). For PP2(k+1), the smallest removal of PP2(k) by bounty
hunters is required coupled with the largest growth of the predator population.
Thus, the correct interval for (24) is:
PP1(k+1) = PP1(k) + DT(PNGR1(k) - PP1(k)4RF2(k)) (30)
PP2(k+1) = PP2(k) + DT(PNGR2(k) - PP2(k)*RF1(k)) . (31)
There are two approaches to converting a system dynamic model into an interval
model when it is represented as in this example, by a list of interconnected
equations. The first way would be to analyze the set of equations carefully to
determine how the endpoints of the interval variables are to be formed as was
done for (30) and (31) above. The second approach would be to state the model
in its vector function representation which for the discrete-time model is:
Keke) = FOR 3K). (32)
Then transform (32) into its interval representation recalling that interval
function evaluation must be performed. That is, (32) becomes:
Xk+1) = F(R(k) 3k) « (33)
The transformation of (8)-(25) into an interval model is given below where the
minima and maxima are taken with respect to the (constrained) state variables
DP(k) © CDP1(k),DP2(k)], PP(k) © CPP1(k),PP2(k)], and F(k) € CF1(k),F2(k)].
DP1(k+1) = min{DP(k) + DT[DNGR(k) - DPR(k)1]}
= min{DP(k) + DT[DP(k)«DGRFT(FR(k)) - PP(k)xDKRT(DD(k))]} , (34)
where FR(k) = F(k)xAFPD/DP(k) and DD(k) = DP(k)/AREA. When average food per
deer AFPD and/or area are themselves positive width intervals, then the minimum
in (34) must also include these two (constrained) variables. Moreover, if or
when the left end of the deer population DP1(k) goes to zero, a suitable
interpretation of FR(k) must be made or the model halted. The right end of the
deer population thus becomes:
DP2(k+1) = max{DP(k) + DT[DP(k)*DGRFT(FR(k)) - PP(k)*DKRT(DD(k))]} . (35)
Likewise,
PP1(k+1) = min{PP(k) + DT[PP(k)*PGRFT(DKR(k)) - PP(k)*RF(k)]} (36)
PP2(k+1) = max{PP(k) + DTEPP(k)xPGRFT(DKR(k)) - PP(k)*RF(k)]} . (37)
When RF(k) e CRF1(k),RF2(k)] is an interval of positive width, then the min/max
above must include it in addition to those of the three state variables. Here
DKR(k) = DKRT(DD(k)) and PGRFT(DKR(k)) is a function of a function; i.e., a
composite function.
F1(k+1) = min{F(k) + DTEGR(k) - DP(k)*FCPDT(FR(k)) J} (38)
F2(k+1) = max{F(k) + DTEGR(k) - DP(k)*FCPDT(FR(k))]} , (39)
wou
-51l1-
where GR(k) = (FCAP - F(k))/FRTT(F(k)/FCAP). As before, when FCAP belongs to
a positive width interval, then the min/max must incorporate these variations.
In (34)-(39) DT, the time increment, is considered as a zero width interval.
What happens then when we transform (8)-(25) "blindly" into an interval model
as illustrated by (27)-(29)? When a straight substitution is made, it would be
equivalent to moving the min/max of (34)-(39) inside the parentheses; that is,
we are obtaining an interval, while guaranteeing to contain all possible
trajectories due to the variations, which is nevertheless far too wide, rendering
potentially "meaningless" results. The interval model (34)-(39) takes min/max
in tandem. For example, consider what occurred in (26) and (28) and compare this
to (36). The minimum in (36) is taken as PP(k) varies over [PP1(k),PP2(k)] (and
the other state variables of course) so that the worse case interval
PNGR1(k) ~ PP2(k)*RF2(k) = PP1(k)«PGRFT(DKR(k)) - PP2(k)«RF2(k) of (28) would
not occur for (36) since PP(k) will be the same value throughout (36). Clearly,
(34)-(39) require a range function analysis or an analytic min/max solution.
However, before computing, the following simplification can/should be made:
DP(k+1) = min{DP(k) + DTLDP(k)*DGRFT(F1(k)*AFPD1/DP(k))
~ PP2(k)aDKRT(DP(k)/AREA2)1} , (34")
where DP(k) ¢ [DP1(k),DP2(k)].
DP2(k+1) = max{DP(k) + DTEDP(k)*DGRFT(F2(k)*AFPD2/DP(k))
- PP1(k)«DKRT(DP(k)/AREA1)]} , (35')
where DP(k) € [DP1(k),DP2(k)].
PP1(k+1) = min{PP(k) + DTLPP(k)xPGRFT(DKR1(k)) - PP(k)«RF2(k)]
= min{PP(k)[1 + DT(PGRFT(DKR1(k)) - RF2(k)) J}
= PP1(k)L1 + DT(PGRFT(DKRI(k)) - RF2(k)I} (36')
where PP(k) © [PP1(k),PP2(k)].
PPp2(k+1) = PP2(k)[1 + DT(PGRFT(DKR2(k)) - RF1(k)]} . (37!)
Equations (38) and (39) do not simplify due to the nonlinearity of GR(k) and
DK(k )aFCPDT(FR(k)).
CONCLUSIONS
An analysis of Tables 1, 2, and 3 indicates that interval analysis is able to
capture the full effect of parameter uncertainty in one pass which would require
the full range of parameter permutations for sensitivity. Secondly, once a set
of equations is put into an interval analytic setting, which in view of (26)-
(29) is not always a one-to-one correspondence between algebraic operations, the
full range of sensitivity analysis can be performed either in the usual fashion
or in a way which analyzes effects due to multiple variations. One is assured
of obtaining the full range of permuted variations in a single analysis.
Secondly, it is clear from the second example, (8)-(37'), that care must be taken
in transforming a model into an interval analytic setting. Whether the effort
involved in doing this is worth the gain accrued in being able to model explicitly
parameter, numerical, and roundoff errors will, of course, depend on the purpose
for which the model serves. However, it has been demonstrated that the methods
of interval analysis can be used to perform just such a task.
REFERENCES
Alefeld, Gétz and Jiirgen Herzberger. Introduction to Interval Computations.
New York: Academic Press, 1983.
-512-
Ford, Andrew, Jeffrey Amlin, and George Backus. "A Practical Approach to
Sensitivity Testing of System Dynamic Models," Proceedings of the 1983
International System Dynamics Conference, 1983, Plenary Session, pp. 261-261.
Forrester, Jay W. and Peter M. Senge. "Tests for Building Confidence in System
Dynamic Models," TIMS Studies in the Management Sciences, Volume 14, 1980,
pp. 201-228.
Forrester, Nathan B. "Eigenvalue Analysis of Dominant Feedback Loops,"
Proceedings of the 1983 International System Dynamics Conference, 1983, Plenary
Session, pp. 178-202.
Goldberg, Samuel. Introduction to Difference Equations. New York: John Wiley
& Sons, 1958.
Goodman, M. R. Study Notes on System Dynamics. Cambridge: M.I.T. Press, 1980.
Graham, Alan K. and Alexander L. Pugh. “Behavior Analysis Software for Large
Dynamo Models," Proceedings of the 1983 International System Dynamics Conference,
1983, pp. 130-142.
Moore, R. E. Methods and Applications of Interval Analysis. Philadelphia:
SIAM, 1979.
Ratschek, H. and J. Rokne. Computer Methods for the Range of Functions. New
York: John Wiley & Sons, 1984,
Tank-Nielsen, Carsten. "Sensitivity Analysis in System Dynamics." In J. Randers
(Ed.) Elements of the System Dynamics Method. Cambridge: M.I.T. Press, 1980.
CC BY-NC-SA 4.0