[From Bill Powers (970121.0100 MST)]
Martin Taylor 970120 11:10 --
One point: my revised manual also notes the computational artifact (found
by Bill P. and me some years ago) involved in the calculation of "optimum
gain" for the "Simple Control System" program, and suggests as an exercise
for the reader that they rewrite the program to simulate the same analogue
network with values of the compute iteration cycle (dt) reduced by a
factor of ten. This augments the discussion of computational artifacts
later in the manual.
That's a good suggestion. I should point out that the hazard of
computational artifacts (which I first discussed in the appendix of B:CP,
and later with more competence in the "Spadework" article in 1979) arose
originally because of the limited speed of digital computers and the problem
of handling huge volumes of output data. When, for practical reasons, you
want to use as few iterations of a program as possible, you make dt as large
as possible, and that's when you run into this problem. Now, with today's
fast computers, this is not so much of a consideration. Wolfgang Zocher has
included, in the latest version of Simcon, the ability to save data only for
every nth iteration (an obvious solution to the data volume problem which I
should have thought of years ago), so there is little remaining reason from
the standpoint of simulation speed or data storage to use large values of
dt. It's a good idea, as you say, to try different values of dt to make sure
that shortening it makes no material difference in the way the program runs,
but there's no longer any need to flirt with values that could create
computational artifacts.
There's a second reason for using small values of dt. The integrations that
are always involved in models of real control systems suffer an
approximation defect, in that the "real" value of the integral is changing
continuously during an iteration, while on the computer it can change only
at discrete intervals. If we were dealing with open-loop systems, the simple
means of integration ("Euler integration") that we use would be completely
inadequate unless we used _extremely_ small values of dt. We would have to
use Runge-Kutta or even more advanced methods of numerical integration in
order to keep the simulated integral from wandering far from the correct
value over even rather short runs of a model (This problem is greatly
reduced in models of closed-loop control systems, because integration errors
do not accumulate; they are handled just like disturbances and are corrected
automatically and invisibly -- if the errors aren't too great).
One way to convince yourself of this is to program a sine-wave oscillator:
x := 100;
y := 0;
{iteration loop}
dx := 0.1*y
dy := -0.1*x
y := y + dy*dt {Euler integrations}
x := x + dx*dt
{end loop}
If you plot x against y, you should get a circle. Making dt as large as the
"optimum" value or even 1/10 of the optimum value will produce a pretty poor
circle -- it will grow or shrink appreciably on every iteration (I forget
which). You have to make dt VERY small (like 0.00001*k) in order for the
circle to repeat itself for more than a few revolutions without any visible
change in size.
The Nyquist sampling criterion says that to preserve all the information in
a sine-wave, you should theoretically have to sample it no more than twice
per cycle. But that has nothing to do with the problem here: here the
problem is to generate a sine-wave that continues for some time and actually
repeats. In general, the problem is to compute the integral of a varying
quantity in a way that will continue to represent the value of the integral
correctly over a long run of the simulation. The value of dt needed to do
this, even in a closed-loop system, is far smaller than the value that would
lead to computational artifacts, and several orders of magnitude smaller
than what the Nyquist sampling criterion would suggest. This is the
principal reason that, in tracking experiments, we have to sample the handle
and cursor positions at least 25 times per second and preferably 60 times
per second or more, even though the bandwidth of a human tracker is between
only 0 to about 2.5 cycles per second. The reason has nothing to do with
"sampling rate." It has everything to do with accurate integration using the
simple method of Euler integration.
In simulating the environment, this is an important consideration. Suppose
we have a position control system, and the object being positioned is a pure
mass on which the system exerts position-independent forces (as in a rocket
ship). The velocity of the mass is the integral of the applied force, and
the position is the integral of the velocity. Any integration errors in
computing velocity are integrated again in computing position, so the errors
are (roughly) squared. To get an accurate model of the environment, you must
either use microscopic values of dt, or use some advanced
predictor-corrector method of integration that minimizes the errors. If you
just use Euler integration with a dt calculated on the basis of the Nyquist
criterion and the maximum expected frequency of variation of position, you
will end up with a horribly inaccurate model of the environment -- too
inaccurate even if the model is a closed-loop control system.
···
--------------------------
Incidentally, there is a related problem I have seen in Hans' way of
discrete modeling. The value of dt is a way of introducing physical time
into a computer simulation. If the model is set up correctly to handle
physical phenomena, then changing dt should change only the simulated
time-scale on which the model runs, not the actual physical time or the
relationships among physical quantities. In effect, the plots look exactly
the same if the same physical time scale is used, with only the density of
plotted points changing as dt is changed (Simcon works this way). In
particular, the model of the plant should run at the same speed in physical
time regardless of the size of dt (as long as it is comfortably shorter than
the minimum value for accurate computation).
However, when the plant is modeled simply in terms of successive states
indexed by k, the speed at which it runs is proportional to the time assumed
for one iteration. The value of dt is always 1 but the units are
unspecified. If you change the association of one iteration with physical
time (in your head), the plant simply runs at a different speed. When you
use indexes that have no specified meaning in physical time, there is no way
to anchor the behavior of simulated physical processes to their actual
physical properties.
Let me use a simple example. Suppose you're simulating the way a mass
converts a force into a velocity. If you simply say
v := v + 0.01*f, in MKS units,
then on each iteration, the velocity will increase by the amount 0.01*f. But
how long a time is represented by one iteration of the loop containing this
expression? If you think of one iteration as representing 0.001 seconds of
physical time, then a force of one newton will create a velocity change of
0.00001 meter per second in one iteration. But if one iteration is imagined
to correspond to 1 second of physical time, then the same force will create
a velocity change of 0.01 meter per second in one iteration.
Obviously, simply thinking of one iteration as representing different
lengths of time CHANGES THE ASSUMED MASS OF THE OBJECT! There is nothing in
the innocent line of program code to suggest this change -- if one weren't
thinking about the physics of what is being simulated, the hidden assumption
of a value of dt, with the subtle effect of changing the modeled mass, could
completely escape notice. In fact, as written and with the units involve, dt
has an implicit value of 1 second, and the coefficient of 0.01 represents
the effect of a mass of 100 kilograms. So you have no choice about the
correspondence between iteration time and physical time: whether you realize
it or not you have assumed a specific value of dt, exactly 1 second.
The correct way to write the program line is
v := v + 0.01*f*dt
where dt is the explicitly assumed duration of one iteration of the program.
Now, if you change dt, you will find that the same physical force always
produces the correct change of velocity in one iteration, assuming a
constant mass of 100 kilograms. For shorter dt's, the velocity changes
appropriately less in one iteration.
Notice that dt enters ONLY in integration steps or other steps including
time functions. It does not enter any intervening calculations involving
only algebraic relationships. This means that some steps of the program
yield different values depending on the assignment of a value to dt, but
other steps do not. That, in turn, means that you can't simply assign
arbitrary meanings to the physical time occupied by one iteration. You must
actually change some of the computations -- but not all -- when you change
the assumed duration of one iteration. The time scale is NOT arbitrary. If
you don't include explicit dt's in the correct places in the program, the
simulation will be correct only for one physical time scale.
I assume that when the successive-state types of models are applied to real
physical situations, the value of dt somehow gets assigned correctly. But I
have yet to see a computation of this kind in which dt appears explicitly.
The general assumption seems to be that the same model, without dt, will
apply to processes that run at any speed.
Best,
Bill P.