Final (?) form of corrected MCT model

[Bill_Powers1]

[From Bill Powers (970306.2143 MST)]

This has been a long and interesting day. Earlier, I posted a preliminary
version of Hans Blom's MCT model, using the correct physical model of the
theodolite. The remainder of the model, however, was uncertain: it
oscillated, and I didn't know why. So I went back to scratch and re-did the
entire model. I am finally satisfied (my, that feels risky!) that the new
version is exactly correct in all particulars, and that I can explain
exactly why it behaves as it does.

For those who aren't going to follow this post in detail, I will summarize
the findings.

1. Except for the matter of the modeled disturbance, the MCT model developed
using Hans' strategy proves identical to the two-level PCT model I used --
and indeed, used in the Little Man some years ago. This was a great
surprise, because the mathematics does not reveal this fact without some
persuasion.

2. For the particular example I chose, in which there is a free-moving mass
accelerated and decelerated in angle by an applied torque, there is an
uncorrectable problem with the MCT model in its present form. Briefly, it is
that the modeled disturbance can't be calculated until after the real
disturbance has imparted an angular velocity to the mass. There is an
unavoidable one-iteration mismatch between the real and modeled
disturbances, and this means that there is an uncorrected angular velocity
in the real theodolite that does not appear in the modeled theodolite. This
mismatch persists as long as there is a nonzero disturbance, and each
attempt to correct it results in another mismatch, so the angular velocity
error is never corrected.

This problem would not be obvious in a physical system that contains
internal damping, as all of Hans' previous models of "the plant" have done,
because the velocity error would quickly die away. I chose this free-moving
mass as an example because I still suspected that slight integration errors
would be cumulative. That proved to be true, although not at all in the way
I had anticipated.

This did not show up in Hans' orginal model, because Hans used identically
the same _incorrect_ calculation for the physical theodolite and the model's
internal simulation of the theodolite. This way of treating the system left
velocity out of consideration, and this prevented the problem I have found.
As it has turned out, the critical error boils down to an uncorrectable
effect of that one-iteration delay between the modeled disturbance and the
real one: truth, it seems, is in the smallest details.

After my haste-induced mathematical "insights" of a couple of days ago, I
don't expect these brave words to be taken without a grain of salt. So I
will spend the rest of this post laying out my development of the MCT model,
after which I hope that more than one person will subject the results to as
careful a scrutiny as I have tried to apply today.

···

---------------------------------------------------------

I: THE THEODOLITE MODEL

Let us start with the model of the theodolite, under the condition that the
torques applied to it are constant except for steps that occur at intervals
of dt, defined for the duration of this development as 0.01 second.

I will briefly recap what I said in an earlier post today.

The basic differential equation of motion is

J*d2x/dt^2 = (u + d),

where x, u, and d are all functions of time. u is the torque output of the
discrete controller, and d is the value of a disturbance that has only two
values: either 0, or some constant number. J is the moment of inertia of the
theodolite.

The solution of this equation in the special case where u and d are constants is

(1) x = x0 + v0*dt + 0.5*(u+d)/J*dt^2

In this equation, dt is not an infinitesimal, but a finite duration of time.
The angle of the theodolite at the beginning of the interval is x0, and the
velocity is v0. The value of x in equation 1 is the position at the end of
the interval dt, and thus at the beginning of the next interval.

Also, we can compute the angular velocity of the theodolite at the end of
the interval. It will be

(2) v = v0 + (u+d)/J*dt

These two equations exactly describe the continuous motion of the theodolite
under the stated conditions. Under torque inputs which occur in discrete
steps, the motion will be a series of segments of parabolas.

II: THE MCT MODEL

The conditions specified by Hans are that the MCT model is to contain a
simulation of the theodolite that is an exact replica of the equations of
motion of the theodolite. This was accomplished by using equations (1) and
(2) for both the theodolite and the simulation.

I chose to start the iteration loop with the torque output initialized to
zero, and compute both the modeled xmod and vmod and the real x and v (using
the term "modeled" to mean "simulated"). The variables a and amod are dummy
variables used to represent angular acceleration. Since these are program
steps, the variable on the left is the _next_ value; the same variable if it
appears on the right is the _old_ value.

Just before each value of x or v is computed, the value of d or dmod is
found. When a disturbance appears, it will appear at the start of an
interval dt.

The first step is to calculate the new position and velocity of the real and
simulated theodolite. We apply the current output u simultaneously to the
simulation and the real theodolite:

   d := make_disturbance; {returns 0 or 50}
   a := (u+d)/J;
   x := x + v*dt + 0.5*a*dt*dt; {real theodolite}
   v := v + a*dt;

amod := (u+dmod)/J;
xmod := xmod + vmod*dt + 0.5*amod*dt*dt; {simulated theodolite}
vmod := vmod + amod*dt;

In the program, the previous values of x and v are preserved as oldx, oldv,
oldxmod, and oldvmod.

Note that in order to calculate the new values for the simulation, we must
already have an assumed value of dmod. However, in order to estimate a new
value for dmod, we must have available the values of x, v, xmod, and vmod,
the latter two of which will be known only after the second set of
calculations above is done. Thus we must do both sets of calculations before
we can make an estimate of the next simulated disturbance.

Since we know that the disturbance begins at zero, jumping to 50 only later,
we can use an initial assumed value of 0.

When the disturbance jumps to 50, we will find that it has been applied to
the real theodolite, but not yet to the simulated one, which is still using
the previous value of dmod. This is where the unavoidable one-iteration lag
of the simulated disturbance arises. The real theodolite receives a brief
impulse that the modeled one does not. As the result of this, the real
theodolite acquires an angular velocity that, as it turns out, remains
higher than that of the modeled system until the end of the disturbance.
This produces a steadily increasing departure of the real theodolite angle
from the simulated angle, for the duration of the disturbance.

After both calculations above have been done, we can estimate the next value
of the disturbance on the basis of any discrepancy between the position and
velocity of the real and simulated theodolites. It is done as follows:

We now have values for x, v, xmod, and vmod all at the same instant of time.
We can thus use the two equations for x and xmod to solve for the mismatch
between the real disturbance and the simulated one. That mismatch can then
be used to correct the simulated disturbance.

   x := x + v*dt + 0.5*a*dt*dt; {real theodolite}

xmod := xmod + vmod*dt + 0.5*amod*dt*dt; {simulated theodolite}

We now have to be careful to keep track of what is a "new" value or an "old"
value. Since x and xmod are new values, obtained from the above two sets of
calculations, we must make use of the old values as they existed before the
calculations, which have been preserved as oldx, oldxmod, and so on. Using
those variables, we can write the above two program steps as equations:

   x = oldx + oldv*dt + 0.5*a*dt*dt (a has not yet changed)

xmod = oldxmod + oldvmod*dt + 0.5*amod*dt*dt

Subtracting the second equation from the first, we have

(x - xmod) = (oldx - oldxmod) + (oldv - oldvmod)*dt
                + 0.5*(a - amod)*dt*dt

Substituting from the definitions of a and amod, we find that the right-hand
term reduces to

0.5*(d - dmod)*dt*dt/J

In that term we find the discrepancy we want to calculate, d - dmod, which
is the amount by which the previous value of dmod must be changed to be a
correct model of the real disturbance.

We can now solve the above equation for (d - dmod):

d - dmod := [(x - xmod) - (oldx - oldxmod) - (oldv -
oldvmod)*dt]*2*J/(dt*dt)

Adding (d - dmod) to the old value, dmod, we get the correct simulated value
of dmod for the next iteration: d. It is one iteration too late, but that is
the best that can be done.

The final step is to calculate the next value of the output, u, preliminary
to starting a new iteration. To do this, we start with the simulation
program step, but using the current values of xmod and vmod in order to
calculate the _next_ value of xmod, which we want to equal the reference
level r.

   xmod(next) = xmod + vmod*dt + 0.5*a*dt*dt

Substituting r for xmod on the left allows us to solve for the next value of
amod that will make the model's xmod equal to the reference value r.

   r = xmod + vmod*dt + 0.5*amod*dt*dt

and expanding amod,

   r = xmod + vmod*dt + 0.5*(u + dmod)*dt*dt/J

Remember that dmod is now the estimate for the _next_ value of simulated
disturbance, already calculated.

From this we obtain the program step for calculating the next value

of u:

u := 2*J/(dt*dt)*(r - xmod - vmod*dt) - dmod

Since we now have new values for u and dmod, we can go back to the beginning
and start the calculations for the next iteration.

----------------------------------------
III: THE PCT MODEL

To clear up one mystery, we must review the PCT model. This is a two-level
model, in which the lower level controls velocity and the higher level
controls position. The error signal from the position control system is
multiplied by Gx to produce Rv, the reference level for velocity. The error
signal from the velocity control system is multiplied by Gv to create the
output torque u. The two control-system equations are

Rv :=Gx*(r - x);
u := Gv*(Rv - v);

The x and v values come from the physical theodolite, which is identical to
the model used for the MCT model above. We assume that the control systems
can sense v and x, sampling once per dt.

The MCT model calculates u this way:

u := 2*J/(dt*dt)*(r - xmod - vmod*dt) - dmod

However, when that equation is used in the MCT model, a serious malfunction
occurs. The output oscillates on alternate iterations between 80 and 0 as
long as the reference signal is changing. The reason is that the error gets
completely corrected on one iteration, leaving no error to drive the next
output, so the error jumps back up again. The cure is quite simple: omit the
2. In fact, the problem here is the same one I spoke of in the "spadeword"
paper: the gain of this system is too high; in fact it is just twice the
maximum gain for correcting the error in minimum time. Reducing the gain by
a factor of two brings the system essentially to its maximum useful
performance (the rapid oscillations disappear). So the program step actually
used in the MCT model is

u := J/(dt*dt)*(r - xmod - vmod*dt) - dmod

The way the next value of u is computed in the MCT model seems to have
little to do with the way it is computed in the PCT model. However, let's
combine the two PCT program steps into one:

u := Gv*[Gx*(r-x) - v)], or

u := Gv*Gx*(r-x) - Gv*v

Let's also multiply out the expression for the MCT model (except for dmod):

u := 2*J/(dt*dt)*(r - xmod) - 2*J/dt*vmod

From this we can compute the relationship between the coefficients in the

two expressions:

    Gv*Gx = J/(dt*dt)

    Gv := J/dt

From this we can deduce that

    Gx := 1/dt

J is 1.0 n-m^2, and dt is 0.01, so the values of Gx and Gv required for the
PCT system to behave like the MCT system without a disturbance are

Gx = 100 and

Gv := 100.

When these values are used in the PCT model, it behaves identically to the
MCT model, including small oscillations when the reference signal suddenly
starts or stops rising.

We're comparing the MCT model controlling its internal simulation with the
PCT model controlling the actual theodolite (since it has no internal
simulation of the theodolite). When there is no disturbance, the two models
behave identically at all rates of rise of the reference signal, all the way
to the point where the slope is so high that both models show large damped
oscillations at the start of the reference signal changes (they could be
overlaid).

So the MCT model controls its _simulated_ xmod in the same way that the PCT
model controls the actual theodolite. But how do the models compare in terms
of controlling the real theodolite?
--------------------------------

With a disturbance of zero, the MCT model controls the real theodolite
exactly as the PCT model does. When a step-disturbance is introduced,
however, the MCT model fails. In the program, the disturbance is introduced
after the reference signal ramp has reached its upper limit and has become
constant. Up to the point where the disturbance begins, the traces for the
PCT model exactly match those of the MCT model. But when the disturbance
appears, the angle of the real theodolite in the MCT model starts to
increase at a constant rate that is proportional to the magnitude of the
disturbance, departing rapidly from the reference signal. This increasing
error is a perfectly linear ramp; nothing is being done to correct it. At
the end of the disturbance square-wave, the increasing angle of the real
theodolite becomes constant again, but does not return to the reference level.

The reason has been explained above: because of the one-iteration delay in
the simulated disturbance, which is unavoidable, the real theodolite gets a
one-iteration kick that the simulated theodolite does not receive. This
causes a velocity to appear in the real theodolite that does not appear in
the simulation. Given the model as it stands, this can't be cured.

The only way to cure it would be to allow the controller to observe the
angular velocity of the real theodolite as well as the angle. Two criteria
would then have to be met in estimating the next disturbance and the next
value of u, and probably the saved previous values of the variables would
have to go back two steps instead of one. Not having tried this, I don't
know what the final model would look like.
---------------------------------------------------
I hope that I have made no errors in this presentation -- I have spent about
16 hours today debugging this analysis, so at least mathematical errors
should have been eliminated. The fact that the PCT model exactly matches the
MCT model in the absence of disturbances strongly indicates that no errors
remain, and the behavior of the theodolite under the MCT model, when a
disturbance appears, seems completely explained.

I hope this presentation will lay to rest some hasty generalizations about
the PCT model. It is neither better nor worse than the MCT model in the main
respects examined here. It is, in fact, the same model, only the PCT model
controls the real theodolite while the MCT model controls an internal
simulation having the same properties as the theodolite. I am sure that a
proper elaboration on the MCT model will take care of the problem of
uncontrolled velocities, and that this will make the two models even more
similar. If this is done, we can then pass on to further items on the agenda:

1. The effects of computational inaccuracies on the two models;

2. The effects of unmodeled disturbances on the two models;

3. Finding a true measure of the complexity of the two models (which, to tip
my hand, must include ALL mathematical operations employed in making the
model work).

I'll post the source code for the final version of theo6mct.pas in the
morning. I can't do any more tonight (or, I guess, this morning).

Best,

Bill P.

[HANS_BLOM2]

------------------- ASCII.ASC follows --------------------
[Hans Blom, 970310f]

(Bill Powers (970306.2143 MST))

1. Except for the matter of the modeled disturbance, the MCT
model developed using Hans' strategy proves identical to the
two-level PCT model I used -- and indeed, used in the Little Man
some years ago. This was a great surprise, because the
mathematics does not reveal this fact without some persuasion.

Well, maybe not "identical", but there is a great "family
resemblance". The reason is, of course, that the best possible
PCT model uses as much knowledge about reality as the MCT model
which has a perfect model.

2. For the particular example I chose, in which there is a
free-moving mass accelerated and decelerated in angle by an
applied torque, there is an uncorrectable problem with the MCT
model in its present form. Briefly, it is that the modeled
disturbance can't be calculated until after the real disturbance
has imparted an angular velocity to the mass.

Yes, this is an unavoidable problem in _any_ model. Something
unexpected and unpredictable happens. Such unforeseen matters
cannot be anticipated, and therefore control cannot be perfect;
it will necessarily lag behind when something unanticipated
occurs. What _can_ be done, however, is to notice what happened
(i.e. the control error) and to correct for it afterwards, so
that it does not have long term effects. This you did not do
correctly:

This mismatch persists as long as there is a nonzero
disturbance, and each attempt to correct it results in another
mismatch, so the angular velocity error is never corrected.

See my original program code where the correction was done
correctly and where there are no cumulative effects.

This problem would not be obvious in a physical system that
contains internal damping, as all of Hans' previous models of
"the plant" have done, because the velocity error would quickly
die away.

You're chasing the wrong rabbit here...

This did not show up in Hans' orginal model, because Hans used
identically the same _incorrect_ calculation for the physical
theodolite and the model's internal simulation of the
theodolite.

Another wrong rabbit. The "perfect" MCT controller expects a
perfect match between model and reality. If you change the way
reality operates (which you do, compared to what my model is
based on), the model is no longer perfect.

This way of treating the system left velocity out of
consideration, and this prevented the problem I have found.

Might it not be true, as an alternative explanation, that you
introduced the velocity -- or whatever -- in an incorrect way?
See, once again, my original program which did not have these
problems, even though it did not perceive velocity.

As it has turned out, the critical error boils down to an
uncorrectable effect of that one-iteration delay between the
modeled disturbance and the real one: truth, it seems, is in the
smallest details.

Errors in control cannot be prevented. They are due to what you
call "disturbances" and what I call "unmodeled dynamics", i.e.
incompleteness of the model. What _can_ be prevented, if we are
lucky, is long range effects of these "disturbances". In my
controller, this was done correctly. In your modification of it,
not anymore, alas.

After my haste-induced mathematical "insights" of a couple of
days ago, I don't expect these brave words to be taken without a
grain of salt.

You're not finished yet, Bill...

I: THE THEODOLITE MODEL

The basic differential equation of motion is

J*d2x/dt^2 = (u + d)

Note that velocity does not appear in the problem formulation. It
is redundant; x can be solved from this equation without
introducing such an "intervening variable": just integrate twice!
(Velocity is the result of integrating once; integrate once more
and it's gone!). Finding the solution (in x) is easy to do for
this differential equation, slightly more involved for the
equivalent difference equation form.

The solution of this equation in the special case where u and d
are constants is

(1) x = x0 + v0*dt + 0.5*(u+d)/J*dt^2

Look for a solution that does not contain a velocity and you're
all set...

II: THE MCT MODEL

This produces a steadily increasing departure of the real
theodolite angle from the simulated angle, for the duration of
the disturbance.

Note that my program didn't have this problem. That suggests that
your re-implementation/reverse engineering is not equivalent ;-).

III: THE PCT MODEL

To clear up one mystery, we must review the PCT model. This is a
two-level model, in which the lower level controls velocity and
the higher level controls position.

To a control engineer this is a strange way to state the problem.
The original specification of the problem was that only position
was to be controlled. Physically it is, of course, impossible to
control both position and velocity simultaneously. If you specify
what one must be, the other is fully determined (except for a
constant of integration). So which aspect(s) of "reality" does
your controller "really" control?

It is no problem, however, to _solve_ the problem by introducing
an extra _computational level_. That's fine with me, except when
this extra level requires an extra (physically and mathematically
redundant!) perception.

We assume that the control systems can sense v and x, sampling
once per dt.

This new assumption was not agreed on, and I find it very
unlogical. It's redundant, and a controller cannot control both
velocity AND position independently. Sensing x has to go...

The cure [of the "MCT" problem] is quite simple: omit the 2.

How did you arrive at this "solution"? By analysis? By trial and
error? Are you sure it is correct? Under all circumstances?

When these values are used in the PCT model, it behaves
identically to the MCT model, including small oscillations when
the reference signal suddenly starts or stops rising. ... When
there is no disturbance, the two models behave identically at
all rates of rise of the reference signal, all the way to the
point where the slope is so high that both models show large
damped oscillations at the start of the reference signal changes
(they could be overlaid).

So far so good...

When a step-disturbance is introduced, however, the MCT model
fails.

It did not in my original code. So it appears that this
observation is a pointer to the place in the "MCT" program code
where something is wrong.

The reason has been explained above: because of the one-
iteration delay in the simulated disturbance, which is
unavoidable, the real theodolite gets a one-iteration kick that
the simulated theodolite does not receive. This causes a
velocity to appear in the real theodolite that does not appear
in the simulation. Given the model as it stands, this can't be
cured.

Bad model, I would say. Correct it...

The only way to cure it would be to allow the controller to
observe the angular velocity of the real theodolite as well as
the angle.

Not "the only way", Bill. Your physical insight ought to tell you
that velocity is completely redundant. Differentiate position and
you have velocity. Integrate velocity and you have position. This
part of the program appears to need more analysis...

I hope that I have made no errors in this presentation -- I have
spent about 16 hours today debugging this analysis, so at least
mathematical errors should have been eliminated.

Not yet, it appears...

The fact that the PCT model exactly matches the MCT model in the
absence of disturbances strongly indicates that no errors
remain

... in the parts of the program that are disturbance-independent.

... and the behavior of the theodolite under the MCT model,
when a disturbance appears, seems completely explained.

May I once more refer you back to my original program that did
not have this problem?

I hope this presentation will lay to rest some hasty
generalizations about the PCT model. It is neither better nor
worse than the MCT model in the main respects examined here. It
is, in fact, the same model ...

One would hope so! Note that you have introduced perfect
knowledge of reality into your controller, just like the MCT
controller has. You have also avoided sub-optimalities such as
"slowing factors" and the like. One would therefore expect that
the controllers behave the same way.

Let's conclude with this observation. I quote "the MCT model
developed using Hans' strategy proves identical to the two-level
PCT model I used -- and indeed, used in the Little Man some years
ago."

That is a fruitful conclusion: we're not so far apart as you
thought previously. So MCT = PCT, after all ;-). And this:

This was a great surprise, because the mathematics does not
reveal this fact without some persuasion.

Right, the mathematics are complex and require analysis of the
finest details. Observe, however, how simple the resulting MCT
program code is. I wouldn't call it significantly more complex
than the PCT code. Would you?

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970310.0811 MST)]

Hans Blom, 970310f --

This mismatch persists as long as there is a nonzero
disturbance, and each attempt to correct it results in another
mismatch, so the angular velocity error is never corrected.

See my original program code where the correction was done
correctly and where there are no cumulative effects.

In your original program code, you used straight-line segments to
approximate the behavior of the theodolite, and the identical model for the
internal simulation in the MCT model. Thus both your model of reality and
the internal simulation misrepresented the actual dynamics of the
environment, and this is why you were able to get "perfect" control using
position information only. When I pointed this out and said we had to use a
phsyically correct model of the theodolite, you rather grudgingly agreed
that this could be done, but said that the internal stimulation still had to
be an exact duplicate of the environment model. I did that.

Another wrong rabbit. The "perfect" MCT controller expects a
perfect match between model and reality. If you change the way
reality operates (which you do, compared to what my model is
based on), the model is no longer perfect.

But the MCT model's internal simulation, in my version, is still a perfect
representation of the environment -- except that now the correct model of
the environment is used. You made the environment fit the MCT model's
incorrect simulation; I made the simulation match a correct model of the
environment. The model is still perfect.

This way of treating the system left velocity out of
consideration, and this prevented the problem I have found.

Might it not be true, as an alternative explanation, that you
introduced the velocity -- or whatever -- in an incorrect way?
See, once again, my original program which did not have these
problems, even though it did not perceive velocity.

It is you who are treating velocity incorrectly. The reason is evident in a
later part of your post:

The basic differential equation of motion is

J*d2x/dt^2 = (u + d)

Note that velocity does not appear in the problem formulation. It
is redundant; x can be solved from this equation without
introducing such an "intervening variable": just integrate twice!
(Velocity is the result of integrating once; integrate once more
and it's gone!). Finding the solution (in x) is easy to do for
this differential equation, slightly more involved for the
equivalent difference equation form.

You will doubtless kick yourself when you recall your basic calculus and see
why this statement is incorrect.

When you integrate twice, you generate two constants of integration. On the
first integration, you have (with v = dx/dt and a = (u+d)/J)

       t = t2
v = integral (a*dt), or
       t = t1

v = v(t1) + a*t

The constant of integration is v(t1).

Integrating again, we get

x := x(t1) + v(t1)*(t2 - t1) + 0.5*a*(t2 - t1)^2

The second constant of integration is x(t1). Note that in this second
integration, the first constant of integration is still present and is
multiplied by t2 - t1.

Physically, this says that the initial position is x(t1) and the initial
velocity is v(t1). The distance travelled due to the initial velocity is
v(t1)*(t2-t1), and the added distance due to the acceleration is
0.5*a*(t2-t1)^2. If the acceleration goes to zero, the theodolite will
continue to swing at the angular velocity v(t1).

The "dt" in the Pascal program is actually (t2 - t1).

The velocity does not "disappear" -- if it does, you have simply made the
mistake of forgetting about constants of integration. When you leave out
velocity, you are claiming that the theodolite moves from the position at t1
to the position at t2 AND STOPS. In fact, it has acquired an increment of
angular velocity which continues during the next interval, adding to the
effects of acceleration during the next interval. This is what you left out
of your model of the theodolite, and of course you left it out of the MCT
model, too.

I think you need to consider this for a while before we continue. I don't
want to go on refuting statements of yours that fail to recognize this
error, and put you on the defensive.

Just to look ahead a bit, if I did calculate velocity in the PCT model
instead of sensing it, the PCT model would then be a functional duplicate of
yours, and THEY WOULD BOTH FAIL TO DEAL WITH THE DISTURBANCE CORRECTLY, for
exactly the same reason. But we can put that subject (an others) off until
you see what I am talking about above. Note that when I say the disturbance
isn't handled correctly, I'm talking about what you see in the lower trace
generated by the program: when the disturbance appears, the position begins
to depart linearly from the reference position at a rather steep angle. This
is not just a brief blip, as in your original program.

Best,

Bill P.

[HANS_BLOM2]

[Hans Blom, 970311b]

(Bill Powers (970310.0811 MST))

In your original program code, you used straight-line segments to
approximate the behavior of the theodolite, and the identical model
for the internal simulation in the MCT model.

I do nothing of the kind. I compute the system's behavior only in a
discrete number of _points_, and I take care that the number of
points is adequate to capture everything "interesting" in the
behavior of the controller. Even the final result, the output on the
computer's screen, is given as discrete _points_. What happens _in
between_ those points is not represented in the difference equation
formalism, is not computed by the control algorithm in the computer,
and is not displayed on the screen. The values _at_ those discrete
time instants, however, are correct, provided only that dt is small
enough.

When I pointed this out and said we had to use a phsyically correct
model of the theodolite, you rather grudgingly agreed that this
could be done, but said that the internal stimulation still had to
be an exact duplicate of the environment model.

Let's recapitulate. The original statement of the problem was

J*d^2x(t)/dt^2 = u(t) + d(t)

The only approximation or transformation that I made was to replace
the second derivative into a "discrete" second derivative, so that we
can work out a digital computer implementation:

    x(t+dt) - 2*x(t) + x(t-dt)
J * -------------------------- = u(t) + d(t)
              dt*dt

Assuming that y'all were not familiar with this transformation, I
provided its derivation in an earlier post. The above formula is what
I take to be the (computational) equivalent of the original second
order differential equation (objections?). The only requirement is
that dt is to be taken small enough. If that is the case, the traces
of an analog computer simulation of the differential equation problem
statement and a digital computer simulation of the difference
equation problem statement could -- for any choices of u(t) and d(t)
-- be overlaid without showing observable discrepancies. In other
words, the choice of dt does not matter as long as it is small enough
(one can even check that using a digital computer only, by varying dt
and overlaying the traces). That is when one could call a model "good
enough", I think.

All else in my MCT controller logically and deductively follows from
this formula ("model") alone.

You made the environment fit the MCT model's incorrect simulation;
I made the simulation match a correct model of the environment.

Please tell me why and where the above transformation is inadequate.

It is you who are treating velocity incorrectly.

Let's see. By rearranging terms, we can transform the difference
equation problem statement

    x(t+dt) - 2*x(t) + x(t-dt)
J * -------------------------- = u(t) + d(t)
              dt*dt

into

x(t+dt) - 2*x(t) + x(t-dt) = [u(t) + d(t)] * dt*dt / J

and thus into

x(t+dt) = 2*x(t) - x(t-dt) + [u(t) + d(t)] * dt*dt / J

and into

x(t+dt) = x(t) + {x(t) - x(t-dt)} + [u(t) + d(t)] * dt*dt / J

and thus again into

x(t+dt) = x(t) + v(t)*dt + [u(t) + d(t)] * dt*dt / J

where v(t) = [x(t) - x(t-dt)] / dt

which is remarkably like your formula. The only difference is that I
compute v from two sequential observations of x rather than measure
it directly. What's wrong?

Greetings,

Hans

[HANS_BLOM2]

[Hans Blom, 970311c]

(Bill Powers (970310.0811 MST))

One more comment:

You will doubtless kick yourself when you recall your basic calculus
and see why this statement is incorrect. ... When you integrate
twice, you generate two constants of integration.

Bill, my controller works with the difference equation directly. The
MCT controller does not integrate. Check (using my original program
code!) what happens when the theodolite starts out with an initial
non-zero velocity. You can do that e.g. by setting xold := -0.1 and x
:= 0.0 in the initialization. You will observe that -- due to its
control actions -- it is rapidly stopped to zero velocity and in the
correct spot. As it should be: it's _position control_. Any deviation
of the velocity will generate a deviation in the position, which will
be controlled away. The same happens with a disturbance: it causes a
velocity change and thus a position deviation, but this deviation
will be controlled away as well.

My original program had no errors, Bill!

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970311.0757 MST)]

Hans Blom, 970311b --

In your original program code, you used straight-line segments to
approximate the behavior of the theodolite, and the identical model
for the internal simulation in the MCT model.

I do nothing of the kind. I compute the system's behavior only in a
discrete number of _points_, and I take care that the number of
points is adequate to capture everything "interesting" in the
behavior of the controller.

Perhaps you didn't _intend_ to do this, but when you use discrete equations
you are always using straight-line segments to replace curves. You then
compute the values at discrete intervals, but the assumption is that between
these intervals the variables are changing continuously and linearly.

Let's recapitulate. The original statement of the problem was

J*d^2x(t)/dt^2 = u(t) + d(t)

The only approximation or transformation that I made was to replace
the second derivative into a "discrete" second derivative, so that we
can work out a digital computer implementation:

   x(t+dt) - 2*x(t) + x(t-dt)
J * -------------------------- = u(t) + d(t)
             dt*dt

You made the environment fit the MCT model's incorrect simulation;
I made the simulation match a correct model of the environment.
It is you who are treating velocity incorrectly.

Let's see. By rearranging terms, we can transform the difference
equation problem statement

   x(t+dt) - 2*x(t) + x(t-dt)
J * -------------------------- = u(t) + d(t)
             dt*dt

into

x(t+dt) - 2*x(t) + x(t-dt) = [u(t) + d(t)] * dt*dt / J

and thus into

x(t+dt) = 2*x(t) - x(t-dt) + [u(t) + d(t)] * dt*dt / J

and into

NOTE: EQUATION NUMBER ADDED

(1) >x(t+dt) = x(t) + {x(t) - x(t-dt)} + [u(t) + d(t)] * dt*dt / J

and thus again into

x(t+dt) = x(t) + v(t)*dt + [u(t) + d(t)] * dt*dt / J

where v(t) = [x(t) - x(t-dt)] / dt

Fine. I have to admit not seeing that your equation does contain a velocity
term -- it's well-concealed, and you did keep saying that we didn't have to
consider velocity. However, this velocity term isn't quite right, because it
is actually v(t - 0.5*dt) -- the average velocity beween t-dt and t.

In that interval, the velocity is continually increasing from v(t-1) to
v(t), so (x(t) - x(t-1))/dt underestimates the velocity at time t.

This underestimation is accomplied by an _overestimation_ of the effect of
acceleration. Comparing your expression with the exact expression, you will
see that a factor of 1/2 is missing from the dt^2 term in your expression.
The effect of acceleration in your expression is overestimated by a factor of 2.

What I want is for us to agree that the expression

x(t+dt) = x(t) + v(t)*dt + 0.5*a(t)*dt*dt

is the _exact_ expression for the motion of the theodolite, given that the
driving force occurs in steps and remains constant during each step. The
appended program compares your expression (1) with mine (2), and with the
exact expression, which is x := 1/2*a*t^2 -- the angle moved in total time t
at an acceleration of _a_. Here _a_ is 10 units, dt is 0.01, and every 10th
calculation is recorded in the output file. xm = mct calculation, xp = pct
calculation, and xe = exact calculation.

t = 0.0000 xm = 0.0000 xp = 0.0000 xe = 0.0000
t = 0.1000 xm = 0.0550 xp = 0.0500 xe = 0.0500
t = 0.2000 xm = 0.2100 xp = 0.2000 xe = 0.2000
t = 0.3000 xm = 0.4650 xp = 0.4500 xe = 0.4500
t = 0.4000 xm = 0.8200 xp = 0.8000 xe = 0.8000
t = 0.5000 xm = 1.2750 xp = 1.2500 xe = 1.2500
t = 0.6000 xm = 1.8300 xp = 1.8000 xe = 1.8000
t = 0.7000 xm = 2.4850 xp = 2.4500 xe = 2.4500
t = 0.8000 xm = 3.2400 xp = 3.2000 xe = 3.2000
t = 0.9000 xm = 4.0950 xp = 4.0500 xe = 4.0500
t = 1.0000 xm = 5.0500 xp = 5.0000 xe = 5.0000
t = 1.1000 xm = 6.1050 xp = 6.0500 xe = 6.0500
t = 1.2000 xm = 7.2600 xp = 7.2000 xe = 7.2000
t = 1.3000 xm = 8.5150 xp = 8.4500 xe = 8.4500
t = 1.4000 xm = 9.8700 xp = 9.8000 xe = 9.8000
t = 1.5000 xm = 11.3250 xp = 11.2500 xe = 11.2500
t = 1.6000 xm = 12.8800 xp = 12.8000 xe = 12.8000
t = 1.7000 xm = 14.5350 xp = 14.4500 xe = 14.4500
t = 1.8000 xm = 16.2900 xp = 16.2000 xe = 16.2000
t = 1.9000 xm = 18.1450 xp = 18.0500 xe = 18.0500
t = 2.0000 xm = 20.1000 xp = 20.0000 xe = 20.0000
t = 2.1000 xm = 22.1550 xp = 22.0500 xe = 22.0500
t = 2.2000 xm = 24.3100 xp = 24.2000 xe = 24.2000
t = 2.3000 xm = 26.5650 xp = 26.4500 xe = 26.4500
t = 2.4000 xm = 28.9200 xp = 28.8000 xe = 28.8000
t = 2.5000 xm = 31.3750 xp = 31.2500 xe = 31.2500

As you can see, in the first step the MCT calculation is off by 10%. The
proportional error gets less as time goes on, but the damage is done in the
first iteration. Note that your expression gives a result that is larger
than mine -- and the difference grows linearly with time at 0.005 per
iteration. Note also that my way of computing x(t+dt) gives the exact value
(and this is independent of the size of dt).

Can you now agree that my way of computing x(t+dt) gives the exact value?

Best,

Bill P.

[HANS_BLOM2]

------------------- ASCII.ASC follows --------------------
[Hans Blom, 970312]

(Bill Powers (970311.0757 MST))

Hans Blom, 970311b --

In your original program code, you used straight-line segments to
approximate the behavior of the theodolite, and the identical
model for the internal simulation in the MCT model.

I do nothing of the kind. I compute the system's behavior only in a
discrete number of _points_, and I take care that the number of
points is adequate to capture everything "interesting" in the
behavior of the controller.

Perhaps you didn't _intend_ to do this, but when you use discrete
equations you are always using straight-line segments to replace
curves.

This just isn't true: _nothing_ is or can be said about what happens
_between_ points; the difference equation formalism only allows us to
talk about values _in_ points.

You then compute the values at discrete intervals, but the
assumption is that between these intervals the variables are
changing continuously and linearly.

There is no such assumption _in the formalism_. Nothing prevents
_you_, however, or anyone else, to entertain such an assumption. The
assumption that the variables change _linearly_ is certainly unsound;
just look at the behavior of the theodolite's position x. In the open
loop case (constant a), the values measured at the discrete sample
times show a quadratic behavior. _If_ you want to entertain some
assumption about how the x-value changes _between_ points, you'd
better use a quadratic interpolation!

If you do want to assume something about what happens between points,
the most useful (and intuitively obvious) assumption is that the
controller's output remains constant in between those points. The
output value is changed, after all, only at discrete time instants,
and that would best correspond with a physical "reality" where the
output remains at the last value it was provided with between those
time instants (but even that depends on the physics of the actuator).
In my MCT theodolite controller, for instance, the output value that
is applied "during" a single period dt and that sets the theodolite
in motion is

u = 200 at dt = 0.001
u = 20 at dt = 0.01
u = 2 at dt = 0.1

where you see that u*dt (the area under the u-curve, i.e. the impulse
given to the theodolite) remains constant and independent of dt. It
is indeed the physically applied impulse that sets the physical
theodolite in motion. We can also easily visualize what would need to
happen if u were limited to 100 at most: at dt = 0.001, the output u
= 100 would need to be given during _two_ consecutive intervals of
duration 0.001 each.

Given the -- almost natural, but not always correct -- assumption
that the controller's output is constant between points (when it is
open loop or "ballistic", after all, and cannot be corrected),
however, it becomes easy to derive what happens between points. This
is, as Martin stated, like _observing_ with a fine comb (a small dt)
while yet _controlling_ with a larger DT >> dt.

My posted theodolite equation could also be written as:

x(t+dt) = x(t) + v(t)*dt + [u(t) + d(t)] * dt*dt / J

where v(t) = [x(t) - x(t-dt)] / dt

Fine. I have to admit not seeing that your equation does contain a
velocity term -- it's well-concealed, and you did keep saying that
we didn't have to consider velocity.

I find "concealed" a pretty subjective notion ;-). But we _don't have
to_ consider velocity. It is present in the formula

  x(t+dt) = x(t) + v(t)*dt + [u(t) + d(t)] * dt*dt / J

but absent in

  x(t+dt) = 2*x(t) - x(t-dt) + [u(t) + d(t)] * dt*dt / J

Yet these formulas _mean_ exactly the same thing, if v(t) is defined
as (x(t) - x(t-dt))/dt. If, however, one defines v(t) as (x(t-dt) -
x(t-2*dt))/dt [another, but slightly worse, approximation], they
don't mean the same thing and will give different numerical results.

However, this velocity term isn't quite right, because it is
actually v(t - 0.5*dt) -- the average velocity beween t-dt and t.

Regrettably, v(t-0.5*dt) is not available; only v(t) and v(t-dt) are.
Unless I would sample at the double rate, of course. In the open-loop
case (a=constant), v(t) gives a slightly too high velocity prediction
and v(t-dt) a slightly too low velocity prediction.

But are you sure that this formula of yours is _always_ correct? Or
just in this specific example, where the acceleration is constant and
the velocity increases linearly? Would the 0.5 need to be different,
maybe, if the velocity increases quadratically? For a quick check,
run the open-loop case with a(t) = t, for instance, where you also
know the analytical solution. I would prefer an analytic derivation,
however.

The difference equation formulation that I used assumed that we only
had measured values available from times t and t-dt, not t-0.5*dt. If
dt is small enough, however, the error does not matter (especially
not in a controller). In your last example, for instance, you
calculate the open-loop solution of d2x/dt2=u, with u constant. My
approximation at time t=2, where xtrue=1/2u*t^2, gives

dt = 0.1 x = 21.000 error = 1
dt = 0.01 x = 20.100 error = 0.1
dt = 0.001 x = 20.010 error = 0.01

An ever better approximation as dt becomes smaller.

In that interval, the velocity is continually increasing from v(t-1)
to v(t), so (x(t) - x(t-1))/dt underestimates the velocity at time
t.

How would you establish the velocity from _measurements_ of x?

This underestimation is accomplied by an _overestimation_ of the
effect of acceleration. Comparing your expression with the exact
expression, you will see that a factor of 1/2 is missing from the
dt^2 term in your expression.

Is this also true in the general case and not just in this example?

What I want is for us to agree that the expression

x(t+dt) = x(t) + v(t)*dt + 0.5*a(t)*dt*dt

is the _exact_ expression for the motion of the theodolite, given
that the driving force occurs in steps and remains constant during
each step.

Before I accept it, can you give a derivation that proves that this
is the correct formula _in the general case_?

Let's go back to the very start. The original formula was (J=1):

  d2x/dt2 = a

We also know that dv/dt = a, so we expect a discrete approximation of
the form

  v(t+dt) := v(t) + a(t)*dt

Would you agree to that, or would you want to modify it to

  v(t+dt) := v(t) + a(t)*dt + ...

where the ... gives an additional correction?

Now the second step. We know that dx/dt = v, so we expect a discrete
approximation of the form

  x(t+dt) := x(t) + v(t)*dt

Would you agree to that, or would you want to modify it to

  x(t+dt) := x(t) + v(t)*dt + ...

where the ... gives an additional correction? You obviously want to
modify it, because you add a correction term 0.5*a(t)*dt*dt. Where
does this term come from? Is it correct _generally_? Why correct the
expression for x but not that for v?

The appended program compares your expression (1) with mine (2), and
with the exact expression, which is x := 1/2*a*t^2 -- the angle
moved in total time t at an acceleration of _a_.

You convincingly demonstrate that my formula is only an approximation
and that your formula is absolutely correct _in this case_: constant
acceleration. Can you also demonstrate that your formula is exact in
the general case? If so, I will have learned something very useful
that will let me improve my models!

I'm sure we can agree on _some_ difference equations of the "true"
theodolite. Given a, we need to have the computer construct x (and,
if you want, v as an intermediate step). You propose the formula

  x(t+dt) = x(t) + v(t)*dt + 0.5*a(t)*dt*dt

In it we need v(t), I presume according to

  v(t) = v(t-dt) + a(t-dt)*dt

Is that correct?

The one fruitful thing that may come out of this whole discussion on
what I think are minor points is that a model-based controller's
model need not be _exact_ for it to yet be pretty good ;-).

The next question that we need to agree on is which perceptions are
available. Is it just position, as I propose, or is it velocity as
well, as you propose? Can you convince me that it is necessary to
measure v? Or is it -- as it appears -- the assumption that it is
impossible to design a _PCT_ controller if v is not available?

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970312.1736 MST)]

Hans Blom, 970312 --

Let's drop the "straight-line segment" argument; it's unimportant. Let's
just skip directly to the main problem here:

Let's go back to the very start. The original formula was (J=1):

d2x/dt2 = a

We also know that dv/dt = a, so we expect a discrete approximation of
the form

v(t+dt) := v(t) + a(t)*dt

Would you agree to that, or would you want to modify it to

v(t+dt) := v(t) + a(t)*dt + ...

where the ... gives an additional correction?

For this special case I will accept the first version. Of course for the
case where u is continuously variable and the relationships are not linear,
we would have to include more terms.

Now the second step. We know that dx/dt = v, so we expect a discrete
approximation of the form

x(t+dt) := x(t) + v(t)*dt

Would you agree to that, or would you want to modify it to

x(t+dt) := x(t) + v(t)*dt + ...

where the ... gives an additional correction? You obviously want to
modify it, because you add a correction term 0.5*a(t)*dt*dt. Where
does this term come from? Is it correct _generally_? Why correct the
expression for x but not that for v?

You are right that I want to modify it. In the first equation, u is in fact
a constant because the controller output is constant between steps.
Therefore the velocity will increase linearly with time, and v2 := v1 +
a*dt, as you have written it. This is exact for the example we're using.

But during this interval, v changes from v1 to v1+a*dt. The distance
traveled is the _average_ velocity (not the final velocity) times dt. The
_average_ velocity is (v1+v2)/2, or

[v1 + (v1 + a*dt)]/2,

yielding

dx/dt := v1 + a*dt/2,

and

x2 := x1 + v1*dt + 0.5*a*dt^2

This is wierd, Hans. All I am doing is going through the elementary calculus
derivation of the double integral of the initial expression. Has it been so
long since you worked with continuous systems that you have truly forgotten
all this? Perhaps now that you see the derivation, it will all come back to you.

You convincingly demonstrate that my formula is only an approximation
and that your formula is absolutely correct _in this case_: constant
acceleration. Can you also demonstrate that your formula is exact in
the general case? If so, I will have learned something very useful
that will let me improve my models!

No, we're not trying to deal with the general case here, which is MUCH more
complex. If you happen to have a solution for the general case, fine -- it
surely should work for this special case. But if it uses an expression for
the theodolite behavior in this special case other than

x2 := x1 + v1*dt + 0.5*a*dt^2

it is wrong, because this is the exactly correct formula for the special
case we are considering.

The next question that we need to agree on is which perceptions are
available. Is it just position, as I propose, or is it velocity as
well, as you propose? Can you convince me that it is necessary to
measure v? Or is it -- as it appears -- the assumption that it is
impossible to design a _PCT_ controller if v is not available?

I tried substituting for the sensed v a very simple calculation of velocity:

v := (x - xold)/dt

which as we know isn't quite correct (one more past value of x would be
needed to get the correct value). The model still worked, although there
were somewhat larger transient oscillations when r and d came to abrupt
change-points. This is to be expected since the actual velocity at time t
will be somewhat larger than computed by this method -- the velocity
feedback is slightly too small and is slightly lagged. The disturbance
caused a slightly-larger departure of x from r while it was present.

We could also calculate the change in v as simply

dv = a*dt,

and add it to a previous v calculated the same way, all the way back to the
start of the run where v is initialized to zero. This would never use a
sensed value of the actual v. Now the calculation of v will be exact, but
for rounding errors.

If we want the exactly correct value of v, and again I repeat _for the
special case we are considering_, we would need x(t), x(t-1), and x(t-2):

[(x(t) - x(t-dt)) - (x(t-dt) - x(t-2dt)]*J/dt = dv/dt at t-1

v at time t would then be v(t-1) + dv/dt*dt.

So this could be done as an independent calculation on each iteration,
rather than being cumulative from the start of the run.

However you calculate the correct value of v(t), you will find that the
calculated v is identical to the v in the theodolite model, and we will be
back to the same PCT program we're already using. So, as it turns out, it is
not critical whether v be sensed or (correctly) calculated: the result will
be the same.

The only difference is that we replace the sensor by a computer that can
receive information about u and J, perform the required calculation, and
output a signal indicating the calculated magnitude of v. I consider that to
be doing it the hard way, when using a tachometer is much cheaper than
adding a computer with the required abilities, considering that the rest can
be done by two 50-cent op-amps.

···

---------------------------------
If I may take it that we agree on the theodolite equation at last, the next
step is to check out the rest of my version of the MCT program.

Since we are using the (new) correct equation to describe the theodolite
movement, the rule is (I presume) that we must use exactly the same equation
as the MCT model's internal simulation of the theodolite. I programmed it
this way:

procedure real_observation;
var savv,savx: real;
{this procedure generates the real x and v}
begin
  savv := v; {Preserving old v and x needed for MCT model only}
  savx := x;
  d := true_disturb(t);
  a := (u + d)/J;
  x := x + v*dt + 0.5*a*dt*dt;
  v := v + a*dt;
  oldv := savv;
  oldx := savx
end;

procedure mod_observation;
var savv,savx: real;
{this procedure generates the simulated x and v}
begin
  savv := vmod;
  savx := xmod;
  amod := (u + dmod)/J;
  xmod := xmod + vmod*dt + 0.5*amod*dt*dt;
  vmod := vmod + amod*dt;
  oldvmod := savv;
  oldxmod := savx;
end;

If you like, you could substitute a calculation of vmod based on past values
of xmod -- as long as they gave exactly the same result.

As you can see, the calculation for the "model observation" (the prediction
of xmod) is exactly the same as that for the "real observation," the actual
position of the theodolite. The prediction of xmod will be exactly correct
at all times if d = dmod. Note that vmod is simply initialized to zero, and
then calculated on each iteration in preparation for the next iteration
where it will be used in the calculation of xmod. The only difference is
that the model observation uses the internally-generated value of
disturbance, dmod, while the theodolite calculation uses the real value, d.

Is this still correct, in terms of the methodology of MCT?

Best,

Bill P.

[HANS_BLOM2]

[Hans Blom, 970318c]

(Bill Powers (970312.1736 MST))

Let's drop the "straight-line segment" argument; it's unimportant.

Is this your modest way of saying you were wrong?

But during this interval [dt], v changes from v1 to v1+a*dt. The
distance traveled is the _average_ velocity (not the final velocity)
times dt. The _average_ velocity is (v1+v2)/2 ...

This is an approximation, which is correct if the velocity increases
linearly but not often otherwise. I can even think of cases where it
is wildly wrong: imagine v going to a very large value right after t1
at which moment we measured v1, yet to come back to the value v2
right before we take the measurement at moment t2. Curing this
problem can only be done by taking smaller values of dt.

This is wierd, Hans. All I am doing is going through the elementary
calculus derivation of the double integral of the initial
expression.

... given some unexpressed assumptions which do not apply in the case
of non-constant acceleration.

Has it been so long since you worked with continuous systems that
you have truly forgotten all this? Perhaps now that you see the
derivation, it will all come back to you.

No, it doesn't. A truly exact "correction" would require an infinity
of extra compensation terms, as you will discover when you try other
circumstances than only a=constant. Try, for instance, the sequence
a=t, a=t^2, a=t^3, etc. By doing so, you would discover the Taylor
series expansion of the correction in terms of dv/dt (a), d2v/dt2,
d3v/dt3 etc. Your correction term is the first of an infinite
sequence.

But if it [a model] uses an expression for the theodolite behavior
in this special case other than

x2 := x1 + v1*dt + 0.5*a*dt^2

it is wrong, because this is the exactly correct formula for the
special case we are considering.

The "special case" we were considering was the theodolite being
controlled. In that "special" case, the acceleration is hardly
constant.

I tried substituting for the sensed v a very simple calculation of
velocity ...

Note that if you want to somehow include your correction term in the
controller, you will need to sense (or compute) the acceleration as
well. Particularly, if the acceleration can change in any possible
way, you will need an infinity of sensors. That is hardly practical.

If I may take it that we agree on the theodolite equation at last

I could agree on almost anything _as a model_ of the real thing. I
cannot agree on your equation as being _exact_, paticularly not if
the theodolite is controlled by a PCT controller ;-).

This strikes me as strange. You _are_ willing to internally model the
acceleration:

  amod := (u + dmod)/J;

but not the velocity. Why is that? Why not sense the acceleration as
well? That would make the problem much easier! ;-).

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970318.1146 MST)]

Hans Blom, 970318--

Let's drop the "straight-line segment" argument; it's unimportant.

Is this your modest way of saying you were wrong?

No. It is a way of saying that I recognize the futility of getting you to
agree on this point, and because it's a peripheral point I think we should
get on to more important matters.

But during this interval [dt], v changes from v1 to v1+a*dt. The
distance traveled is the _average_ velocity (not the final velocity)
times dt. The _average_ velocity is (v1+v2)/2 ...

This is an approximation, which is correct if the velocity increases
linearly but not often otherwise. I can even think of cases where it
is wildly wrong ...

Of course you can, but none of those cases is the one we're considering. I
am fully aware of what happens to the solutions of differential equations
when the driving function is made general: the usual result is that there is
no analytical way to solve the equations. If this were an analog model the
problem would be much more difficult to handle analytically. However, we are
speaking of discrete controllers, which can vary their outputs only in
steps, with the output being constant between steps. This is about the tenth
time that I have carefully spelled out this condition.

This is wierd, Hans. All I am doing is going through the elementary
calculus derivation of the double integral of the initial
expression.

... given some unexpressed assumptions which do not apply in the case
of non-constant acceleration.

These assumptions are NOT unexpressed. I keep expressing them and you keep
ignoring them. The driving torque from the discrete controller IS CONSTANT
BETWEEN STEPS. However, the theodolite, being a real physical device, can
move only continuously. We are calculating the theodolite's continuous
motion under any pattern of step-changes in torque occurring at intervals of DT.

Has it been so long since you worked with continuous systems that
you have truly forgotten all this? Perhaps now that you see the
derivation, it will all come back to you.

No, it doesn't. A truly exact "correction" would require an infinity
of extra compensation terms, as you will discover when you try other
circumstances than only a=constant.

That is totally irrrelevant: we are dealing with a particular case. In the
imaginary case you proposed, where velocity is proportional to acceleration
times time, we have another special case, and solving this one does not
require an infinite number of terms, either. You need an infinite number of
terms only when you try to deal with every possible case, and we are not
doing that here. Your MCT model can't deal with every possible case, either:
there are cases in which there are no inverses available, so you can't solve
the critical equation for the disturbance or the output. We can deal with
that later, if you want to. Right now, we're just trying to handle one
simple special case.

But if it [a model] uses an expression for the theodolite behavior
in this special case other than

x2 := x1 + v1*dt + 0.5*a*dt^2

it is wrong, because this is the exactly correct formula for the
special case we are considering.

The "special case" we were considering was the theodolite being
controlled. In that "special" case, the acceleration is hardly
constant.

The special case so far is only that of a theodolite with no friction or
restoring spring being driven by a torque motor which steps from one torque
to another. During each period between steps, the torque, and therefore the
acceleration, is constant, permitting us to calculate the continuous motion
of the theodolite until the next torque step occurs. At that point we have a
new torque and acceleration, and we can calculate the motion during the next
interval. This calculation is exact whether the duration of DT is 0.01
second or 100 seconds.

I tried substituting for the sensed v a very simple calculation of
velocity ...

Note that if you want to somehow include your correction term in the
controller, you will need to sense (or compute) the acceleration as
well. Particularly, if the acceleration can change in any possible
way, you will need an infinity of sensors. That is hardly practical.

No, the PCT model does not require a sensor for the output or the
acceleration -- although with such a sensor, less well-behaved physical
plants can be handled better.

If I may take it that we agree on the theodolite equation at last

I could agree on almost anything _as a model_ of the real thing. I
cannot agree on your equation as being _exact_, paticularly not if
the theodolite is controlled by a PCT controller ;-).

The model of the theodolite is independent of the model of the controller. I
use exactly the same model of the theodolite in the PCT and the MCT
controllers. We're just modeling a physical object subject to forces.

This strikes me as strange. You _are_ willing to internally model the
acceleration:

amod := (u + dmod)/J;

but not the velocity. Why is that? Why not sense the acceleration as
well? That would make the problem much easier! ;-).

In the MCT model I do model both the acceleration and the velocity. In the
above statement that struck you as strange, I was speaking of substituting a
modeled velocity for the sensed velocity in the _PCT_ model. And I showed
you that if a _correct_ model of the velocity is used, the result is exactly
the same as if the real velocity were sensed. So the issue of sensing is a
non-issue. The issue, if any, is how accurate a calculated model can be.
Since we are both assuming, for the time being, that any internal simulation
of the environment is perfect, this issue is not relevant.

Best,

Bill P.

[HANS_BLOM2]

(Bill Powers (970318.1146 MST))

------------------- ASCII.ASC follows --------------------
[Hans Blom, 970319]

(Bill Powers (970318.1115 MST))

I wish you'd stop calling this a "correction term". It is simply the
right way to represent the solution of the differential equation. I'm
not making it up just to get a lucky right answer; I'm going through the
calculus manipulations.

But let's not loose sight of their physical meaning!

Your new program uses a new acceleration proportional to time. This
requires a new solution of the differential equations, which I have
tediously worked out (not being very good at this); you have to find
the integral and then evaluate it from time t to t + DT, and now we're
dealing with cubic equations. After filling about eight sheets of paper
trying to get it right, I ended up with the appended program, which
again gives the exact position.

I have some problems with your approach. Not that the computation is not
correct (it is, assuming that ...), but that it is an unlikely _model_.
Look at your formulas

  x2 := x2 + v2 * dt + (t*dt*dt + dt*dt*dt/3.0)/2.0; {MY NEW X}
  v2 := v2 + (2.0*t*dt + dt*dt)/2.0; {MY NEW VELOCITY}

and note that this is not the standard difference equation formulation
anymore. For one thing, you assume the controller to have an exact inner
clock that accurately delivers the time t. I have no doubt that organisms
(we want to go there, don't we?) do have internal clocks, but those are
hardly this precise. If the clock isn't exact, these formulas introduce a
severe form of what you have called "integration error". Moreover, these
formulas -- even with an exact clock -- are of no use _in a controller_:
they are fully bound to the postulate that a=t and cannot be generalized
to something that is useful for control purposes.

My "correction term" applies to the real theodolite, and is necessary to
represent its response to the controller's output correctly.

It does not apply to any "real" theodolite, which is far more complex; it
is a (severely simplified) _model_ of a theodolite; see Martin's
comments. But let's skip that for the moment. What I do hope you realize
is that you also assume some very unrealistic properties of the "real"
_controller_. This time, however, these are spelled out explicitly. That
helps.

(Bill Powers (970318.1146 MST))

Here you spell out your control assumption in full detail.

... we are speaking of discrete controllers, which can vary their
outputs only in steps, with the output being constant between steps.

This may have to do with the earlier "straight line" confusion. Note how
unrealistic this assumption is, physically. It requires the actuator to
go from one output level to another one in zero time; in practice, in a
time negligible compared to dt. That would be a reasonable assumption/
approximation if dt is large (which it cannot be for another reason --
there is no active control _during_ dt) but not for the values that are
used in a controller.

This is about the tenth time that I have carefully spelled out this
condition.

You have not. You may have assumed it and taken it for granted, but you
have never spelled out how unrealistic an assumption _about the actuator_
you accept while insisting on utter "realism" (exact matching of
differential and difference equation model) for the theodolite itself.

On the other hand, I have always insisted and taken for granted that we
only have an _approximation_ of the real thing and that such is fine: if
the approximation is good enough, good enough control will result; where
the "good enough" becomes ever better for smaller values of dt (errors in
computations allowing) -- but at the cost of increased computation time.

At the meta-level: we have both interpreted our exchanges in the light of
our existing "internal models", which -- obviously -- severely color what
we "perceive" ;-).

Here is your assumption again:

The driving torque from the discrete controller IS CONSTANT BETWEEN
STEPS. However, the theodolite, being a real physical device, can move
only continuously.

And the controller and its actuator are _not_ real physical devices? It
can _immediately_ go from one torque value to another? I find this
inconsistent. And it is probably this utter inconsistency that has led me
"astray" all this time in not understanding you.

In the imaginary case you proposed, where velocity is proportional to
acceleration times time, we have another special case, and solving this
one does not require an infinite number of terms, either. You need an
infinite number of terms only when you try to deal with every possible
case, and we are not doing that here.

In a "real" controller with a "real" actuator, we _do_ have "every
possible case". But now that I have understood you we need not quarrel
about this anymore, even though we disagree about what meaningful models
are.

Your MCT model can't deal with every possible case, either: there are
cases in which there are no inverses available, so you can't solve the
critical equation for the disturbance or the output. We can deal with
that later, if you want to.

Right. I have always understood your emphasis on "disturbances" and your
insistence that control can be good despite them as an indication that
you might be able to understand that _approximate_ models might result in
good control as well. MCT controllers do not depend on having an exact
model. In practice, an exact model isn't possible anyway. Yet, control
will be good. But we'll come to that later, indeed.

The special case so far is only that of a theodolite with no friction or
restoring spring being driven by a torque motor which steps from one
torque to another. During each period between steps, the torque, and
therefore the acceleration, is constant, permitting us to calculate the
continuous motion of the theodolite until the next torque step occurs.
At that point we have a new torque and acceleration, and we can
calculate the motion during the next interval. This calculation is exact
whether the duration of DT is 0.01 second or 100 seconds.

All right, that's explicit enough. I finally got it...

... the PCT model does not require a sensor for the [velocity?] output
or the acceleration -- although with such a sensor, less well-behaved
physical plants can be handled better.

All right, let's go on to the next question. I am willing to agree with
you on having an unrealistic controller/actuator that can change its
output torque _immediately_ however small dt is. That assumption will,
indeed, allow an external observer to compute what happens _between_
steps, if one wants, even though the controller will not perform such
calculations.

So what do we measure? Position only? That was my assumption right from
the start, and I suggest it even now because I think it is the most
realistic: in practice, the velocity sensor adds too little information
to be worth the expense. Can we agree on that? We might have to
relinquish some "exactness", but that is an illusion anyway.

The model of the theodolite is independent of the model of the
controller.

In what way? My MCT controller presupposed that a "perfect model" was
available, i.e. that the controller knows the equations of motion of the
theodolite, i.e. how values of u and d result in a value of x. You use
those in your controller as well. So we seem to agree on that. Do we?

So what do you mean by the above? I can only understand it in the sense
that in our respective controller implementations the "theodolite data
generator program code" has "hidden variables" that the "controller
program code" cannot access but has to reconstruct if it requires their
(approximate) values [reminds you of The Test?]. In our case -- if we
measure position only, as I suggest -- these "hidden variables" would be
velocity, acceleration and disturbance. The "controller program code"
_is_ allowed to use 1) the form of the theodolite's equations of motion
(I'll accept your version) and 2) the knowledge of the actual theodolite
object itself (in our case that would be only J). Can we agree on that as
well?

The "theodolite data generator program code" is to be based on the
following formulas:

a(t) = [u(t) + d(t)] / J

This one seems obvious. Next, is v(t) to depend on a(t) or on a(t-dt),
i.e. which of the two following formulas for v do we use?

v(t) = v(t-dt) + a(t) * dt
v(t) = v(t-dt) + a(t-dt) * dt

And then, which formula for x do we use? These are the possibilities, I
presume:

x(t) = x(t-dt) + v(t) * dt + 0.5 * a(t) * dt * dt
x(t) = x(t-dt) + v(t) * dt + 0.5 * a(t-dt) * dt * dt
x(t) = x(t-dt) + v(t-dt) * dt + 0.5 * a(t) * dt * dt
x(t) = x(t-dt) + v(t-dt) * dt + 0.5 * a(t-dt) * dt * dt

where u is the controller's output (constant during dt) and d is the
disturbance (zero initially and finally, constant during some episode).

Note that all are good models (i.e. approximations!) for control purposes
if dt is small and changes are not too abrupt. Pick the ones you prefer
based on your considerations of "exactness", and I'll agree with you --
regardless of physical considerations such as infinitely fast effects in
the actuator. Our next task would to be to establish that the program
code correctly implements those formulas. You probably prefer to go the
other way and derive the formulas from the program code ;-). Yes, I know
I'm going slowly...

For the moment, let's assume that the "controller program code" knows J
and dt (both are constant) and u, and observes x (at every iteration) but
not d, a or v. Can we agree on that?

I hope we're finally converging...

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970319.0904 MST)]

Hans Blom, 970319 --

I wish you'd stop calling this a "correction term". It is simply the
right way to represent the solution of the differential equation. I'm
not making it up just to get a lucky right answer; I'm going through the
calculus manipulations.

But let's not lose sight of their physical meaning!

That is exactly what I am trying to avoid.

I have some problems with your approach. Not that the computation is not
correct (it is, assuming that ...), but that it is an unlikely _model_.
Look at your formulas

  x2 := x2 + v2 * dt + (t*dt*dt + dt*dt*dt/3.0)/2.0; {MY NEW X}
  v2 := v2 + (2.0*t*dt + dt*dt)/2.0; {MY NEW VELOCITY}

and note that this is not the standard difference equation formulation
anymore. For one thing, you assume the controller to have an exact inner
clock that accurately delivers the time t.

These formulas have NOTHING T0O DO WITH THE CONTROLLER. They are a physical
model of the assumed theodolite presented with step-torque inputs from ANY
source. If you shrink dt arbitrarily close to zero, they become the
equations of continuous motion as a physicist would derive them. The
controller can sample this behavior at regular or irregular intervals DT;
how the position (and velocity, if used) are sampled is the controller's
problem. If the controller samples at intervals DT, you could run 10 or 100
iterations of the above steps during each interval, with dt = 0.1 or 0.01
times DT, and use DT intervals independent of the dt intervals (to the
nearest whole dt).

My "correction term" applies to the real theodolite, and is necessary to
represent its response to the controller's output correctly.

It does not apply to any "real" theodolite, which is far more complex; it
is a (severely simplified) _model_ of a theodolite; see Martin's
comments.

If you mount a theodolite tube on the rotor of a pancake-style
low-inductance torque motor, the equations of motion will match the above
real behavior within a tenth of a percent or better. I don't know where
you're getting this "severely simplified" idea.

This may have to do with the earlier "straight line" confusion. Note how
unrealistic this assumption is, physically. It requires the actuator to
go from one output level to another one in zero time; in practice, in a
time negligible compared to dt. That would be a reasonable assumption/
approximation if dt is large (which it cannot be for another reason --
there is no active control _during_ dt) but not for the values that are
used in a controller.

That't not true. The torque developed by a torque motor can, with suitable
driving electronics, change from one level to another in 20 microseconds.
The main limitation is how large an inductive voltage spike the driver
transistors can withstand.

This is about the tenth time that I have carefully spelled out this
condition.

You have not. You may have assumed it and taken it for granted, but you
have never spelled out how unrealistic an assumption _about the actuator_
you accept while insisting on utter "realism" (exact matching of
differential and difference equation model) for the theodolite itself.

I certainly have mentioned this condition before, many times. You have
simply not seen its import. Unless you plan for your controller to take its
samples and perform all its computations 50,000 times per second, this
assumption is not at all unrealistic. If you do plan to sample and compute
that fast, then you need to write the model of the theodolite so it takes
into account the smooth inductance-limited change in torque in response to a
step-input of command voltage (which can change in 10 nanoseconds or less).

On the other hand, I have always insisted and taken for granted that we
only have an _approximation_ of the real thing and that such is fine: if
the approximation is good enough, good enough control will result; where
the "good enough" becomes ever better for smaller values of dt (errors in
computations allowing) -- but at the cost of increased computation time.

That depends on where the approximation is. The ultimate test of the
controller would be to have it drive a real theodolite and detect its
position using a signal from a real potentiometer and perhaps a tachometer.
Then the theodolite itself would be doing these "calculations." I am simply
trying to develop a model of the theodolite that would reproduce its actual
behavior as closely as we could measure it, for the particular case we are
considering. This is only a substitute for having a physical theodolite
mounted on a torque motor (with drive electronics) to play with. If we had
the real system we could dispense with modeling it (at least as far as the
PCT model is concerned).

Here is your assumption again:

The driving torque from the discrete controller IS CONSTANT BETWEEN
STEPS. However, the theodolite, being a real physical device, can move
only continuously.

And the controller and its actuator are _not_ real physical devices? It
can _immediately_ go from one torque value to another? I find this
inconsistent. And it is probably this utter inconsistency that has led me
"astray" all this time in not understanding you.

Let us say that the torque motor can go from one torque to another in 20
microseconds. This involves no instant motion of the theodolite; it merely
changes the acceleration that fast, by changing the current through the
stator windings in 20 microseconds. Since we've been assuming values of dt
from 0.001 seconds upward, this means that the torque can change from one
value to another in at most 1/50th of one interval dt. The motions of the
theodolite, however, are smooth and continuous. The only discontinuity is in
the second derivative.

The controller, while being a perfectly physical device, is being
implemented in a digital computer, and can produce only discrete outputs
which are constant between changes. I see no inconsistencies here: both
devices are being analyzed in terms of their physical characteristics.

···

-------------------

All right, let's go on to the next question. I am willing to agree with
you on having an unrealistic controller/actuator that can change its
output torque _immediately_ however small dt is. That assumption will,
indeed, allow an external observer to compute what happens _between_
steps, if one wants, even though the controller will not perform such
calculations.

Fine. I'll let that "unrealistic" pass, in the interests of getting on with
this discussion.

So what do we measure? Position only? That was my assumption right from
the start, and I suggest it even now because I think it is the most
realistic: in practice, the velocity sensor adds too little information
to be worth the expense. Can we agree on that?

I guess you don't get my point that it makes no difference whether you sense
velocity or compute it, as long as we're assuming that the internal
simulation matches the equations of motion exactly. We're not trying to save
a few cents here; the vast cost of sensors and their utter unreliability are
figments of someone's imagination, anyway. The control electronics are far
more likely to fail, just because there are so many more components
(especially if they have to do all the computations you're proposing).

The PCT model uses a velocity sensor as well as a position sensor. If you
don't want to use one, fine. If you do, that's fine with me, too. You design
your model, I'll design mine.

The model of the theodolite is independent of the model of the
controller.

In what way? My MCT controller presupposed that a "perfect model" was
available, i.e. that the controller knows the equations of motion of the
theodolite, i.e. how values of u and d result in a value of x. You use
those in your controller as well. So we seem to agree on that. Do we?

Yes, but that is only to get past this initial hurdle of how to model the
theodolite. In the real MCT model you have to consider modeling errors, too
-- the internal simulation will never be exactly like the real external
world. We will get around to that question, and when we do we will need have
a model of the theodolite that reflects its actual behavior, not its
behavior as simulated inside your controller. We will have to consider, for
example, rounding errors in the calculations.

The "controller program code"
_is_ allowed to use 1) the form of the theodolite's equations of motion
(I'll accept your version) and 2) the knowledge of the actual theodolite
object itself (in our case that would be only J). Can we agree on that as
well?

Yes.

The "theodolite data generator program code" is to be based on the
following formulas:

a(t) = [u(t) + d(t)] / J

Yes.

This one seems obvious. Next, is v(t) to depend on a(t) or on a(t-dt),
i.e. which of the two following formulas for v do we use?

v(t) = v(t-dt) + a(t) * dt
v(t) = v(t-dt) + a(t-dt) * dt

The second, which can also be written

v(t+dt) = v(t) + a(t)*dt.

Just to prevent confusion, let's agree on the convention that t represents
the states of the variables at the _beginning_ of the interval dt.

And then, which formula for x do we use? These are the possibilities, I
presume:

x(t) = x(t-dt) + v(t) * dt + 0.5 * a(t) * dt * dt
x(t) = x(t-dt) + v(t) * dt + 0.5 * a(t-dt) * dt * dt
x(t) = x(t-dt) + v(t-dt) * dt + 0.5 * a(t) * dt * dt
x(t) = x(t-dt) + v(t-dt) * dt + 0.5 * a(t-dt) * dt * dt

The fourth is, I believe, the correct one, but check me on this. At the
beginning of the previous interval (time = t - dt), we have the initial
position x(t-dt), the initial velocity v(t-dt), and the acceleration which
is constant starting at t - dt and remains the same until time t. The
increment in x added to the initial value is the initial velocity times dt,
plus the distance traveled due to the acceleration: 1/2 a*dt^2. Here dt
represents a finite length of time. The same relationship can be expressed as

x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt

After we calculate the next x, we can calculate the next v according to

v(t+dt) = v(t) + a(t)*dt.

The acceleration a, of course, is (u + d)/J.

Our next task would to be to establish that the program
code correctly implements those formulas. You probably prefer to go the
other way and derive the formulas from the program code ;-). Yes, I know
I'm going slowly...

For the moment, let's assume that the "controller program code" knows J
and dt (both are constant) and u, and observes x (at every iteration) but
not d, a or v. Can we agree on that?

As you wish, for the MCT model. I am responsible for the PCT model, and will
propose any sensors I think are necessary. The PCT model does not need
information about J or d, although in effect it has to be manually adjusted
to take J into account (at this point -- later on I will propose a
higher-level system that adjusts parameters for the main system).

I hope we're finally converging...

We'll see. Perhaps at this point you could lay out the basic steps for the
MCT model, after which I will check them against my program code to see if I
got them right (you can double-check me). I'm not sure that I calculated the
next u using information from the right time-step, and ditto for the
prediction of d.

Best,

Bill P.

[HANS_BLOM2]

[Hans Blom, 970320]

Bill, the program below will duplicate your results; for now, the
"controller" outputs a constant level. Check out procedure "observe".
It implements the theodolite's motion equations.

Is this all right?

Greetings,

Hans

program theodolite_model;

const
  J = 1.0; {theodolite characteristic}

var
  x, {perceived position}
  v, {perceived velocity}
  u, {controller's output}
  r, {reference level}
  dt, {time increment; var so it can be changed}
  t: real; {time}
  c: integer; {counter}

procedure observe;
{this procedure implements the theodolite motion equations; it
generates the NEXT observation(s) x(t+dt) and v(t+dt), given u(t)
and d(t); thus it also shifts the time by dt; note that aa, vv and
xx are LOCAL variables, not visible outside this function; declaring
them as "const" saves their value from run to run; aa, vv and xx are
initially known to be zero}

  function disturbance: real;
  {this function implements the disturbance, which is known only
   to procedure observe, the theodolite's motion simulator; this
   function is declared within procedure observe so that it will
   not be visible elsewhere}
  begin
    disturbance := 0.0 {dummy for now}
  end;

const
  aa: real = 0.0;
  vv: real = 0.0;
  xx: real = 0.0;
begin
  aa := (u + disturbance) / J;
  xx := xx + vv*dt + 0.5*aa*dt*dt;
  vv := vv + aa*dt;
  t := t + dt;
  x := xx; v := vv {the controller may use these}
end;

function reference (t: real): real;
{this procedure specifies the reference level r(t) for x(t)}
{not used yet}
begin
  if t < 3.0 then
    reference := 0.0
  else
  if t < 6.0 then
    reference := (t - 3.0) / 3.0
  else
    reference := 1.0
end;

procedure control;
{this function implements the controller; it knows the motion equa-
tions and J, is allowed to observe x(t) and v(t), and delivers u(t)}
begin
  u := 1.0 {dummy for now}
end;

begin
  c := 0; {counter}
  dt := 0.01; {or any other value}
  t := 0.0; {start at t=0}
  x := 0.0; {initially zero}
  v := 0.0; {initially zero}
  repeat
    if (c mod 10) = 0 then {show results}
      writeln (t:12:6, v:12:6, x:12:6);
    r := reference (t+dt); {specify desired x}
    control; {controller computes u}
                                  {time shifts by dt here}
    observe; {observe x (and v?)}
    inc (c)
  until t > 2.001
end.