[From Bill Powers (970301.1450 MST)]
Hans Blom sent me his preliminary MCT model for the theodolite example,
asking me to add the graphics and return it, which I did. The MCT model
oscillated because of the limits put on the output, so I sent it back to
Hans to allow him to find the problem and fix it. However, this discussion
is of interest to CSGnet, so I think that from here on we should conduct the
exchanges in public.
To recap, I will reconstruct Hans' development of his model, mostly to make
it clearer in my own mind. This is a long derivation for which I apologize,
but if there are any conceptual errors they are very subtle and the only way
to show that they do or do not exist is to slog through the details.
Let r(t) = reference signal at time t
x(t) = angle at time t
a(t) = acceleration at time t
v(t) = velocity at time t
u(t) = torque at time t
D(t) = estimated disturbing torque ("disturbance") at time t
(using a cap D to avoid confusing D(t) with dt)
J = moment of inertia of theodolite
dt = time increment
The basic differential equation is
(1) J*d2x/dt^2 = u(t) + D(t)
To compute the second derivative on the basis of observed angles to date, we
take the difference between two velocities (first derivatives). Note that
the following is a little different from what Hans posted, although the
result is the same:
x*t) - x(t - dt)
(2) v(t - 0.5*dt) = -----------------
dt
x(t-dt) - x(t - 2*dt)
(3) v(t - 1.5*dt) = ----------------------
dt
The acceleration at time t-dt is the change in velocity divided by dt:
v(t-0.5*dt) - v(t - 1.5*dt)
a(t - dt) = -----------------------------, or
dt
[x(t) - x(t-dt)] - [x(t-dt) - x(t - 2*dt)]
a(t - dt) = ----------------------------------------------, or
dt
x(t) - 2*x(t - dt) + x(t - 2*dt)
(4) a(t - dt) = ------------------------------------
dt*dt
This acceleration at time t - dt, times the moment of inertia J, is equal to
the total applied torque at time t - dt:
x(t) - 2*x(t - dt) + x(t - 2*dt)
J* ------------------------------------ = u(t - dt) + D(t - dt)
dt*dt
NOTE: at the time this equation holds true, we have observed the present
value of x(t), because t always denotes present time. However, we cannot
write the equation for the acceleration at time t, because it requires
calculations involving previous and future values of x. The acceleration we
calculate is always what existed on the _previous_ iteration, and thus the
equation always refers to the _previous_ values of u and D.
This equation can be used to estimate the value of x(t+dt) given u(t) and
D(t). Hans did this by first adding dt to indices on both sides of the equation:
x(t+dt) - 2*x(t) + x(t-dt)
J* ------------------------------------ = u(t) + D(t).
dt*dt
Solving for x(t+dt) given u(t) and D(t), we obtain the _prediction_
(5) x(t+dt) = [u(t) + D(t)]*dt*dt/J + 2*x(t) - x(t-dt).
The strategy of the MCT model is then to say that x(t+dt) is to be equal to
the reference signal r(t+dt), and to compute the value of u(t) that will
give this result. If x(t+dt) = r(t+dt), we then have
r(t+dt) = [u(t) + D(t)]*dt*dt/J + 2*x(t) - x(t-dt),
which we can solve for u(t). This is the _present_ value of u(t) which will
yield x(t+dt) = r(t+dt).
(6) u(t)= [r(t+dt) - 2*x(t) + x(t-dt)] * dt^2/J - D(t);
NOTE that x(t) is the result from the _previous_ iteration and the
_previous_ value of output, u(t-dt). The value of u(t) remains to be
estimated, and D(t) has not yet been established; when they are established,
they will determine the _next_ value of x, x(t+dt).
The initial value of D(t) is set to zero, and from then on it is found by
looking at the difference between x(t+dt) as predicted, and the value of
x(t+dt) obtained by observing the plant. The key to understanding the method
is that the "observation" of x uses an equation exactly the same as equation
(5), except that the true disturbance is used instead of the estimated
disturbance D(t). Any discrepancy between the predicted and actually
observed values of x must be due to an incorrect value of the
estimated D(t). Thus the assumed value of D(t) can be changed in the right
direction and by the right amount, at most one iteration late, on the basis
of the difference between the observed and predicted values of x(t+dt).
Obviously, if the method of predicting x(t+dt) does not exactly reflect the
input-output function of the plant, the estimated disturbance will be
adjusted automatically to make up for the difference. In the present case
the equations for the plant are identical to those used for the prediction,
so this problem does not arise.
···
-----------------------------------------
However (and I am afraid I must impose on your patience even more, here):
Just to see what would happen, I substituted equation (5) for the method of
integration used to convert torque into angular velocity and position in my
model. It didn't work! The model ran, but when I changed the value of dt,
the basic time increment, the model became slower or faster to respond,
instead of showing exactly the same response plotted with a different
density (as it should). This also happens in Hans' program -- the
oscillations it currently suffers change in frequency when dt changes, and
the time taken to reach the final state changes.
Just as a cross-check, I substituted MY method of computing the theodolite
response to torque for Hans' method, in his program. Now Hans' program
didn't work. In addition to the oscillations caused by limiting the output,
there was a persistent oscillation of output that didn't ever go away.
Perhaps Hans will find some simple mistake that will cure these problems. I
can't see what it would be, since equation (5) is derived in a
straightforward way. But I hope to hear from Hans soon (on CSGnet).
I don't want to post Hans' program until he has had an opportunity to fix
any errors in it, particularly those that cause oscillations when the output
is limited. When there is no limit to the output, the model works perfectly,
correcting the error due to a sudden change of reference signal in one
iteration. However, the torque required to do this, with dt set to 0.001
sec., is 250,000 newton-meters. Considering that the limit of 750
newton-meters is equivalent to the maximum torque of a large truck engine,
the output produced by Hans' model is somewhat ridiculous.
There is much more to say, but I will wait until Hans releases his model for
posting on CSGnet.
Best,
Bill P.