[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.