System Dynamics '95 — Volume I
Chaos from Generic Structures:
A Cautionary Tale
George P. Richardson
Rockefeller College of Public Affairs and Policy
University at Albany — State University of New York
Albany, NY 12222
Abstract
An elegantly simple and understandable model composed of generic structures is shown to move
from damped oscillations to sustained oscillations to repeated period doubling, leading eventually
to chaos. The third-order structure contains a balancing loop with a delay and integral control that
attenuates as the accumulated pressure gets relatively large. The model and its behavior are so
classic and it so quickly converges on its attractor that it may well serve as an ideal structure with
which to introduce complex nonlinear behavior and deterministic chaos. Unfortunately, the model
contains a subtle error, which, when corrected, destroys the model’s ability to exhibit complex
behavior. Reflection on these results suggests a number of cautions about modeling practice.
Introduction
The model described in this note is interesting for two reasons. First, it is an elegant blend of well-
known generic structures that unexpectedly exhibits complex nonlinear dynamics and deterministic
chaos. Furthermore, it does so remarkably quickly: it moves to its various attractors rapidly
enough to make it an excellent choice for introducing complex nonlinear dynamics and
deterministic chaos in real time in class situations or other demonstrations. In light of the model’s
extremely nice features, the second reason for the model’s interest is most unfortunate. The model
contains a subtle but instructive flaw (which the reader is invited to try to find before reading the
last part of this note) which, when cured, kills the model’s chaotic tendencies. The purpose of this
note is to put forward the model both as an elegant and speedy introduction to the phenomena of
period doubling and deterministic chaos, as well as an instructive cautionary tale with perhaps far-
reaching implications.
Model structure
The model diagramed in Figure 1 contains two generic structures, a balancing loop with a delay
(loop 1; Senge 1990) and a second negative loop developing integral control pressure (loop 2).
The balancing loop with delay represents the common effort to reach and maintain a goal in a
context of information delays. The integral control pressure represents the common tendency to
increase efforts to return the system state to its goal the longer the goal-gap persists. Increasing
corrective action the longer a system fails to attain its goal adds considerable instability. Yet
integral control is common: a daily example is our tendency to press the return key repeatedly
when our computer fails to respond quickly enough. More seriously, integral control tendencies
probably characterize many policy areas. Long ago, for example, Phillips (1954) suggested that
monetary policy may have integral control tendencies, as when the Federal Reserve continues to
raise interest rates to cool a stubbornly overheated economy.
The level for accumulated pressure is formulated with an outflow to allow the accumulation to
dissipate over time if the current pressure for change remains small. The integral control term in
the rate of change in the current state is modulated by an effect of extremes of pressure (loops 3
and 4) that brings the weight on accumulated pressure to zero when it grows as large as 50 percent
of the perceived state. (See Figure 1 for a graph of the Effect of Extremes of Pressure.) That
effect was crafted to represent the idea, thought common in political situations, that accumulating
pressure for results eventually comes to be ignored by those who believe they are taking proper
211
Plenary Program
action to solve a problem.!
Model equations are shown in Figure 2. The model is initialized in equilibrium (current state =
perceived state = goal = 100; accumulated pressure = 0) and disturbed at time 3 by a step in the
goal to 120. Simulations that follow were run with a DT of .25 using fourth-order Runge-Kutta.
Model behavior
With a perception time of 6 months and an adjustment time of 3 months in the current pressure for
change, the balancing loop with delay (loop 1) by itself generates damped oscillations (not shown)
with a period of about 25 months. Thus with the normal weight on accumulated pressure set to
zero, the system has a point attractor. As the normal weight on accumulated pressure is increased
from 0 to 0.5, the behavior of the model moves through a series of more complex behaviors and
attractors.
With the normal weight on accumulated pressure set to 0.1, the oscillations are still damped.
When the weight is 0.2, the oscillations expand and become stable; the attractor thus shifts from a
point to an oval orbit. The behavior becomes significantly more interesting for larger values of the
weight.
Figure 3 shows the times series behavior and state-space plots of the model simulated with the
normal weight on accumulated pressure set to 0.35. After a very short transient (about one and a
half cycles in the state-space plots, or 18 months in the time series), the behavior settles into the
high-peak/low peak pattern characteristic of period-2 behavior. As a consequence of the high-
peak/low-peak pattern, the period has doubled to about 50 months.
Figure 4 shows the corresponding plots that result from a normal weight on accumulated pressure
of 0.37. The behavior moves through a transient lasting about half the simulation and then settles
into an apparent period-4 pattern which is more evident in the time series graph than the state-space
plots where the transient tends to obscure the eventual attractor. (The transient period can be
speeded up considerably and the period-4 pattern made more vivid; see the appendix on notes for
real-time demonstrations.) The period has doubled again, to about 100 months.
The plots in Figure 4 do not by themselves guarantee that the simulation has settled on its attractor,
or what its attractor really is. Further information can be obtained by approximating the largest
Lyapunov exponent in the system from the graph shown in Figure 5. It is a defining characteristic
of chaotic models that they exhibit extreme sensitivity to initial conditions. To determine if this
simulation run has such a sensitivity, an exact duplicate of the model can be formulated and run
together with the original, with one level in the duplicated model given a slightly different initial
condition (Wolf 1985; Chen 1988). Here, accumulated pressure was initialized in the second
model at 0.001 instead of 0.0. One then computes the Euclidean distance between the two models?
and plots its natural logarithm over time, as in Figure 5. The slope of that plot is an approximation
of the largest Lyapunov exponent of the system. If it is positive, the system is exhibiting
deterministic chaos. If it is zero, as in Figure 5, the system is periodic or quasi-periodic.
The graph shown in Figure 5 contains a bit more information. It appears to settle into a narrower
range of variation after t = 500 , suggesting that the model has settled onto its attractor after t=
500. To see the attractor more clearly, one can create state-space plots for the later portion of a
1 The effect of extremes of pressure is formulated as a function of the ratio of the absolute value of accumulated
pressure to the perceived state; the absolute value enables the formulation to capture both positive and negative
extremes of pressure.
2 Compute the Euclidean distance as [(x1i— x2i)2 ]¥2, where the sum is taken over all pairs of corresponding
levels x1 and x2 in the two models.
212
System Dynamics '95 — Volume I
run.3 Figure 6 shows such a state-space plot for the simulation in Figures 4 and 5, showing data
for t= 600 to 1200. The plot suggests that the system has indeed settled onto its attractor, but it is
a period-8 attractor (eight orbits in closely arranged pairs), not period-4. When this particular
model is simulated with Euler integration, the transients are much shorter, enabling this sort of
behavior to be seen much more easily; see the appendix.
The plots shown in Figure 7 show that deterministic chaos is evident when the normal weight on
accumulated pressure is set at 0.39. The time series plot (run for 1200 months instead of 800)
shows no repeating pattern, the state-space plots are blurs of nonrepeating orbits, and, more
significant, the graph whose slope approximates the largest Lyapunov exponent has a positive
slope.4 Thus for the weight set to 0.39 we have obtained deterministic chaos.
The result is remarkable in that it is achieved in a simple model consisting of two generic structures
or archetypes: a balancing loop with a delay, coupled with accumulating (integral) control
pressure. The nonlinear effect that attenuates the integral control pressure when it becomes large
relative to the perceived state is apparently crucial to the structure’s ability to exhibit complex
nonlinear dynamics. The complex nonlinear dynamics of this simple structure add impetus to the
focus on dynamic complexity rather than detail complexity in policy modeling (Senge 1990). The
plausibility of the structures in the model suggest that it might be fruitful to try to document
en policy structures having this attenuating integral control in a system with a tendency to
oscillate.
A flaw, and its implications
This rather lovely result is undermined by a flaw in the model which was unnoticed by the author,
and a number of experienced modelers he presented it to, for far too long. Has the reader detected
it?
The flaw stems from an unfortunate choice in the defaults in the simulation software STELLA II
(and iThink). It is startling to realize that if the model equations shown above were simulated in
any other simulation language available to system dynamics modelers today, the results would not
match the results obtained here. And the differences would be major: the equations in any other
language will not produce complex nonlinear dynamics and chaos. Now has the reader located the
problem?
The flaw appears in this model in the outflow from accumulated pressure, formulated to allow that
pressure to dissipate over time when the current state remains near the goal. The unfortunate
default in STELLA II and iThink is to make that outflow a one-way flow — it dissipates
accumulated pressure as long as the pressure is positive, but shuts the flow to zero when the level
goes negative. In this case, negative pressure is meaningful (pulling down the current state of the
system when it is too high for too long), and dissipating negative pressure is also meaningful. The
model, as interpreted by STELLA II, does not dissipate accumulated pressure when it is negative.
Sadly for chaos enthusiasts, when the flow is corrected in STELLA II to be a bi-flow, the ability of
the model to exhibit complex nonlinear dynamics disappears.5
3 In STELLA II one must create dummy variables that are set equal to constants for t < t* and are set equal to the
state variables for t > t*. Scatter plots of these dummy variables will show the the state-space plots from t = t* on.
4 The drop around t = 900 and subsequent behavior saturating under y = 4 are due to the diameter of the attractor:
the points orbiting each attractor in the two versions of the model have become essentially as far apart as possible.
To be certain of the chaotic regime, one can periodically bring the points close together again and continue the
simulation to see the positive slope of [X(x1i— x2j)2 ]12 repeatedly reasserted.
5 I believe it would be a serious mistake to conclude that all deterministic chaos in system dynamics models comes
from flaws in model formulation although Saeed (1992) has argued in that direction. See, e.g., Sterman (1988).
213
Plenary Program
The model thus serves to illuminate several important cautions for system dynamics models, some
cautions about software and some cautions about searching for complex nonlinear dynamics.
Cautions about software
The default in STELLA II that fooled me (and others, I’m anxious to note) for far too long here is a
potential trouble spot for all modelers. One-way flows are intended to prevent levels from going
negative when they should not: population, inventories, bank balances (usually), and so on.
Making outflow rates shut down when levels drop to zero is essential for such levels. The
teaching literature in the field consistently talks about first-order control of outflow rates for levels
that conceptually can not take on negative values (e.g., Richardson and Pugh 1981, 139;
Richmond, Vescuso, and Peterson 1987). The one-way flow in STELLA II is an attempt to
capture in’software an automatic first-order control, so the modeler can ignore it.
Unfortunately, modelers do ignore it (as we see), to their peril. The default one-way flow control
can be activated in a model without the modeler’s knowledge. The behavior of the model can be
influenced by a hidden structure, not made explicit in the conceptualization and formulation
process. This author has seen students fail to catch modeling flaws because they were hidden by
STELLA II’s one-way flow, and he should have known better in his own modeling, but he was
caught anyway. But there is an even more telling caution about STELLA II’s one-way flow
default: leaving the control of outflows to an automatic choice in software allows the modeler to
omit and ignore important control processes that must exist in the reality being modeled. I believe
peas particular software default is an enemy of good conceptualization and good modeling
abits.
Thus, the first caution we learn from this attempt to draw chaos from generic structures is to pay
attention to the software defaults. A reasonable thing for experienced modelers to do using
STELLA II or iThink is to make every flow a bi-flow (or go back to using STELLA 2.10). One
might suggest that it is a reasonable thing to teach to beginning modelers as well, but then one
wonders why we should be bothering with such a counterproductive default in one of the principal
software packages in the field. That observation leads to the next caution.
Cautions about model translation
The model listing above could not be transported to other simulation software packages and
achieve the same observed dynamics. The one-way flow is not perceivable in the equations. (The
translation routines in S**4 and Vensim, for example, would not pick up the one-way flow; the
translated model would not show complex nonlinear dynamics.) To implement such a one-way
flow in other languages, one would need to employ a MAX function, a CLIP, or an IF-THEN
statement, and nothing in the listing alerts us to the need.
Thus a second caution we learn is that translation among software packages may be difficult, and
we are alerted to look for hidden software elements that must be explicitly formulated in other
languages. Given that there are now a plethora of available simulation environments for system
dynamicists — S**4, DYSMAP, Professional DYNAMO, STELLA I/iThink, COSMIC, Vensim,
PowerSim, some of which are available on more than one type of machine or operating system —
this translation problem could become acute. One hesitates to suggest it, but do we need field-wide
minimum standards for software to insure that modeling work can be replicated in diverse
environments?
Cautions about looking for chaos
It is reasonable to suggest that I was drawn into a mistake in this work because of my excitement in
seeing period doubling in a model based on simple generic structures. The promise of ending a
modeling life of unrequited love for complex nonlinear dynamics could be said to have clouded my
perception. I was consumed, one might say, by a lust for chaos.
214
System Dynamics '95 — Volume I
Make no mistake: the lure of chaos is powerful. Not very long ago an otherwise ordinary study
became grist for Science or Nature if the model on which it is based could be drawn into chaotic
regions. I would indeed have loved to present this work as complex nonlinear dynamics from
generic structures, without any modeling flaws. So would others: everyone I present this story to
quickly gravitates to thinking about how the structure can be made meaningful in spite of the one-
way flow on the dissipation of accumulated pressure. Maybe, for example, the accumulated
pressure to reduce criminal justice budgets does not dissipate if they are needlessly high, while the
pressure to push up budgets when they are low does have the dissipative tendency. Such
speculations are theory in search of reality. Although I engaged in them myself in the course of
this work, I believe they have no place in good model-based analyses of real problems.
So I derive another caution from this tale: while one should be alert for complex dynamic
tendencies in models, one should be very careful to ground them in data and to investigate them
thoroughly using realistic parameters. It is no longer interesting to produce an unrealistic model
that can exhibit chaos. We must not be seduced by the siren song of complex dynamics.
Cautions about finding chaos
The original work with this model was done using Euler integration, and the original motivation
for writing about it was the remarkable speed the model settled on its attractors. Although Euler
integration is nowhere near adequate for serious investigation of complex nonlinear dynamic
patterns, I still believe that this model is a superb one to use to introduce the important notions of
attractors, period doubling, sensitivity to initial conditions, estimating the largest Lyapunov
exponent from a graph, and the like. However, the analyses using fourth-order Runge-Kutta
show that previous conclusions were no quite so well founded. Runge-Kutta simulations do not
converge to the system’s attractors as quickly, once the key parameter is greater than 0.35. And,
not surprisingly, the parameter values associated with particular patterns such as period-4 or
period-8 behavior are different in the Runge-Kutta simulations.
The caution we novice chaos explorers learn here is to use the best integration schemes available to
us. But we learn more. This model also shows us that it is hard to nail down a particular complex
dynamic pattern. Figure 4 shows how easy it is to be misled by simulation runs that are too short
or too lacking in resolution to allow a definitive conclusion. We learn we must be extremely
skeptical of the character of dynamic patterns emerging from a model capable of exhibiting
complex nonlinear dynamics.
References
Chen, P. 1988. Empirical and Theoretical Evidence of Economic Chaos. System Dynamics
Review 4(1-2): 81-108.
Phillips, A.W. 1954. Stabilization Policy in a Closed Economy. Economic Journal (June 1954):
290-305. See Richardson (1991, 136-139).
Richardson, G.P. 1991. Feedback Thought in Social Science and Systems Theory. Philadelphia:
University of Pennsylvania Press.
Richardson, G. P. and A. L. Pugh IIL. 1981. Introduction to System Dynamics Modeling with
Dynamo. Portland, Oregon: Productivity Press.
Richmond, B., P. Vescuso, and S. Peterson. 1987. An Academic Users Guide to STELLA.
Hanover, NH: High Performance Systems.
Saeed, K. 1992. Chaos out of Stiff Models. Proceedings of the 1992 International System
Dynamics Conference, Utrecht, The Netherlands, July 1992.
Senge, P.M. 1990. The Fifth Discipline: the Art and Practice of the Learning Organization. New
York: Doubleday/Currency.
Sterman, J. D. 1988. Deterministic Chaos in Models of Human Behavior: Methodological Issues
and Experimental Results. System Dynamics Review 4(1-2): 148-178.
Wolf, A., J. Swift, H. L. Swinney, and J. Vastano. 1985. Determining Lyapunov Exponents
from a Time Series. Pysica D 16: 285-317.
215
Plenary Program
Appendix - Notes on demonstrations
Though too inaccurate to allow reliable conclusions about complex nonlinear dynamics, Euler
integration dramatically improves presentations of the essential dynamics of this model. The model
moves to its attractor usually very quickly and exhibits the same sequence of patterns of behavior
as the Runge-Kutta runs, although for somewhat different values of the key parameter, the normal
weight on accumulated pressure. Given that the more accurate simulations using fourth-order
Runge-Kutta have been performed, it seems sensible to use Euler integration in presentations. As
a further aid to speeding the results, one can change the middle value of the table function for the
multiplier on the weight on accumulated pressure from 0.945 to 1; this tiny change divides the
table into two linear sections and appears to sharpen the model’s movement to its attractors in both
the Runge-Kutta and Euler simulations. The runs in this paper, including those described below,
use the arguably more realistic table containing a curve instead of a corner.
Table 1 is a quick guide to what one can find using Euler integration with this model. Runs were
performed in STELLA 2.2.1 on a Mac IIsi (68030 with 68881 FPU) and a Mac Centris 650
(68040); given the numerical sensitivity of chaotic systems it might be reasonable to expect that
results might differ in other software or on other platforms. Note that in order to obtain complex
nonlinear dynamics with this model in other software one must formulate the dissipating rate from
accumulated pressure to shut down to zero when the level is negative.
Weight Behavior Comment
0.200 Limit cycle | Period = 25 months
0.280 Period 2 Converges very quickly
0.300 Period 2 Converges very quickly
0.320 Period 4 Converges very quickly
0.330 Period 4 Converges very quickly, well separated orbits
0.332 | Period 4 --> 8?|Still fuzzy after t > 600
0.333 Period 8 Transition from period 4 to 8 complete by t = 600
0.334 Period 8 Converges almost immediately
0.335 ? Long transient
0.336 Chaos Lyapunov graph shows three dips, but generally rises
0.340 Chaos After t = 600, looks like period 4, but Lyapunov exponent > 0
0.350 Chaos Lyapunov exponent > 0; slope always > 0
0.400 Chaos Converges immediately to 5 vivid orbits, but drifts to chaos
0.500 Chaos Starts _with 10 orbits
Table 1: Behavior of the model using Euler integration, tabulated against selected values of the
normal weight on accumulated pressure.
216
System Dynamics '95 — Volume I
Figures
Figure 1: Model
diagram. Loop lisa
balancing loop with
a delay in the
perceived state.
Loop 2 is an integral
control loop. Loops
3 and 4 interact to
shut down the Net incr in cum pressure (Loop 1)
integral control effect
if accumulated d ed state
pressure becomesa (3Q===—=== Current state
significant fraction
of the perceived ;
state. The graph Forgetting cum pressyre
below shows the
assumed Effect of
Extremes of
Pressure.
fessure
2 of extremes of pressure
spegnestiit
Accumulated_pressure(t) = Accumulated_pressure(t - dt) + (Net_incr_in_cum pressure -
Forgetting_cum_pressure) * dt
INIT Accumulated_pressure = 0
Net_incr_in_cum_pressure = Current_pressure_for_chng
Forgetting_cum_pressure = Accumulated_pressure/12
Current_state(t) = Current_state(t - dt) + (Chng_in_current_state) * dt
INIT Current_state = Goal
Chng_in_current_state = Current_pressure_for_chng+Wt_on_cum_pressure*Accumulated_pressure
Adj_time = 3
Current_pressure_for_chng = (Goal-Perceived_state) /Adj_time
Extremity_of_accumulated_pressure = abs (Accumulated_pressure) /Perceived_state
Goal = 100 + step(20,3)
Nrmal_wt_on_cum_pressure = .38
Perceived_state = smthl (Current_state, 6)
Wt_on_cum_pressure = Nrmal_wt_on_cum_pressure*Efct_of_extremes_of_pressure
Efct_of_extremes_of_pressure = GRAPH (Extremity_of_accumulated_pressure)
(0.00, 1.00), (0.05, 1.00), (0.1, 1.00), (0.15, 1.00), (0.2, 1.00), (0.25, 1.00), (0.3,
.945), (0.35, 0.79), (0.4, 0.525), (0.45, 0.265), (0.5, 0.00)
Figure 2: Model listing. The model is disturbed from equilibrium by a step in the Goal at t = 3.
217
Plenary Program
Accumulated pressure
30.00
i |
Ly i
ih
:
i
i i
Figure 3: Time series
behavior and state-
space plots of the
model simulated with
the normal weight on
accumulated pressure
set to 0.35. Period-2
behavior.
.00 200.00
180.00
Perceived state
120.00 =
400.00
Months
60.00
0.00 . 125.00 .
Current state
220.00
Current state
80.00
. 4
(
a
3
&
Zz 0.004
3
3
3
2
-80.00 s " 1 i
6000 = * 120.00 "480,00
218
Perceived state
System Dynamics '95 — Volume I
Accumulated pressure
80.00
80.00 my
0,00 =
- i |
FIG negrrneeeresennmomesnsteenseeceits ope initio To
| |
at
Te
et
.00 200.00
Figure 4: Time series
behavior and state-
space plots of the
model simulated with
the normal weight on
accumulated pressure
set to 0.37. Apparent
period-4 behavior in
the time-series plot
after a transient of
about 500 months.
Perceived state
Accumulated pressure
125.00 “220.00
Current state
10.00
219
180.00
120.00 =
125.00
0.00 mf
80.00
120.00
Perceived state
180.0
Sj...
Plenary Program
Figure 5: Graph of the 8.00
natural log of the
Euclidean distance
between two versions
of the model with very
slightly different initial
conditions, for the
simulation shown in
Figure 4. Zero slope
indicates periodic or
quasi-periodic
behavior.
0,00 at
loge({ Z(xsi — x2i)2 1*/2))
a
Figure 6: State-space plots for
the simulation shown in
Figures 4 and 5, showing data
for t = 600 to 1200 to eliminate
the transient patterns. Period-
8 behavior is discernable from
the eight closely paired orbits.
Accumulated pressure
300.00 900.00 1200.06
180.00 =
Perceived state
128.00
Current state
10.00
Accumulated pressure
Current state
Perceived state
System Dynamics '95 — Volume I
B06
125.0031
0.00 300.00 600.00 SORTS 1200.00
Months
80.00
Figure 7: Plots of the model
simulated with the normal
weight on accumulated
pressure set to 0.39, fort =O 7
to 1200. Deterministic chaos
is indicated by the positive
slope of In([X(xui - x2i)2]1/2). a
Accumulated pressure
80.00
In((Zbx1i - x2i)2]1/2)
221
