Improved theodolite program

[Bill_Powers1]

[From Bill Powers (970305.0501 MST)]

Hans Blom (970504...) and other programmers --

Attached is theo5mct.pas, which is a rearrangement of theo4mct.pas. The main
changes involve moving all environmental equations into "do_observation",
which is now a procedure instead of a function.

In "do_observation," it will be seen that there are _two_ environmental
equations, one used for the MCT model and the other for the PCT model. There
is a boolean variable, "ismct", which is now set by whichever model is
running, so the correct environmental model can be chosen. By inserting a
"not" before "ismct," you can verify that these environmental models work
only with the respective controllers. Therein lies a problem for the MCT model.

In the MCT model there is a "predictor" equation which is an exact
derivation from the environmental model. However, the environmental model is
not exactly, physically, correct. When the physically correct environment
model is used with the MCT model, the output of the MCT model oscillates
continuously -- although the controlled variable behaves essentially the
same as before. Control is now achieved by varying the duty cycle of a
square-wave oscillation between maximum positive and negative output at the
sampling frequency.

The physically correct environmental model is a simple application of the
laws of motion. If a constant torque u is applied to the theodolite for a
time dt, the angle x will vary as

x := x0 + v*dt + 0.5*a*dt*dt

The velocity v will increase as

v := v + a*dt.

The angular acceleration a is just (u + d)/J: torque divided by moment of
inertia. Since (u + d) is constant during the interval dt, this is the exact
way to calculate the value of x at the start of the next iteration.

This being the case, both the MCT and the PCT models should use the exact
equation of motion, and not the approximation in the original MCT model. An
adjustment has to be made to the "predictor" equation to correct for the
problems that come from using the exact environmental equations in the MCT
model. I will leave it to Hans to make the adjustment.

Appended is the new program; I will delay posting it with the runnable
version to my FTP page until Hans has the MCT model running with the correct
environmental model.

Best,

Bill P.

program theo5mct;
{
This program compares the performance of the MCT and PCT models. Both
models are controlling the position of a theodolite, under the
following conditions:

change of position: 1 radian.
time resolution: 0.01 sec (dt).
Maximum available output torque: 100 newton-meters
Moment of inertia of theodolite: 1 newton-meter^2
Rate of change of reference signal: variable from keyboard
disturbance: 100 n-m, 3 sec duration, random starting time.

The MCT model runs first, in the lower half of the screen. The slope
of the change in reference signal can be adjusted up and down (by
factors of 1.5) by pressing the + and - keys, shifted or unshifted. A
new run is done every time any key is pressed -- except 'p' or 'q'.

Pressing the 'q' key exits the program at any time. Pressing the 'p'
key switches to running the PCT model; after the first time, pressing
'p' or 'm' switches control to the PCT or MCT model, for further
adjustments of the slopes.

The beginning time of the disturbance varies randomly each time a new
run occurs. By pressing the space bar, new runs can be done with the
same value of the slope of the reference signal's rise but a different
start of the disturbance. Thus you can see what happens when there
is or is not a disturbance occurring during the change in reference
signal.

MCT program by Hans Blom; PCT program and presentation by W. T. Powers
4 March 1997
}

uses
  dos, crt, graph, grutils;

var
  J, K, dt, r, x, u, a, v, maxu, xold, xpre, xsav: real;
  d, d1, d2, trued, t, pslope,mslope: real;
  gv, gx, rv, rx: real;
  xplot: integer;
  maxx, maxy,ycenter: integer;
  ch: char;
  numstr: string;
  fo: text;
  ismct: boolean;

function make_reference (t: real; var slope: real): real;
  {this function defines the setpoint at time t}
var ref: real;
begin
if t >= 3.0 then
begin
  ref := slope*(t - 3.0);
  if ref > 1.0 then ref := 1.0;
  make_reference := ref;
end
  else make_reference := 0.0;
end;

function true_disturb (t: real): real;
  {this function defines the true disturbance}
begin
  if (t < d1) or (t > d2) then
    true_disturb := 0.0
  else
    true_disturb := 50.0;
end;

procedure do_observation;
{this function generates the x that the controller will observe}
begin
  trued := true_disturb(t);
  xsav := x; {save present x}
  if ismct then
   begin
    x := 2.0 * x - xold + u/K + trued/K;
   end
  else
   begin
    x := x + v*dt + 0.5*a*dt*dt;
    v := v + a*dt;
   end;
  xold := xsav; {xold := previous x}
end;

procedure plotit(baseline: integer);
begin
  xplot := round(50.0*t)+ 160;
  putpixel(xplot,baseline,white);
  putpixel(xplot,baseline - round(100.0*r),yellow);
  putpixel(xplot,baseline - round(100.0*x),white);
  putpixel(xplot,baseline - round(trued),lightred);
  putpixel(xplot,baseline - round(u),lightcyan);
end;

procedure legends;
begin
  clearviewport;
  setcolor(white);
  outtextxy(0,0,'Pointing angle');
  setcolor(yellow);
  outtextxy(0,15,'Ref level');
  setcolor(lightcyan);
  outtextxy(0,30,'Output torque');
  setcolor(lightred);
  outtextxy(0,45,'Disturbance');

end;

Procedure PCTmodel;
begin
ismct := false;
setviewport(0,0,maxx,ycenter,ClipOn);
repeat {REPEAT WHOLE RUN}
  t := 0.0; {start at zero time}
  d := 0.0; {assume no disturbance initially}
  v := 0.0;
  x := 0.0;
  gv := 100.0;
  gx := 50.0;
  j := 1.0;

  d1 := 1.5 + 6.0*random;
  d2 := d1 + 3.0;
  maxu := 0.0;
  legends;
  str(pslope:5:3,numstr);
  setcolor(white);
  outtextxy(0,60,'Slope = ' + numstr);
  outtextxy(150,100,'PCT MODEL');

repeat {the control loop starts here}

  r := make_reference (t+dt,pslope); {define reference}

  {COMPUTE OUTPUT}
  rv := gx* (r - x); {velocity ref level = output of position control}
  u := gv*(rv - v); {output force = output of velocity control}

  {LIMIT OUTPUT TO +/- 100 N-M}
  if u < -100.0 then u := -100.0 else {limit output, if desired}
    if u > +100.0 then u := +100.0;

  {ENVIRONMENTAL EQUATIONS}
  do_observation;
  plotit(130);
  t := t + dt;

until t >= 9.0; {at this point the loop ends}

ch := readkey;
if ch in ['=','+'] then pslope := pslope*1.5;
if ch in ['_','-'] then pslope := pslope/1.5;
until ch in ['q','Q','m','M'];
end;

Procedure MCTmodel;
begin
ismct := true;
setviewport(0,ycenter+1,maxx,maxy - 20,ClipOn);
repeat {REPEAT WHOLE RUN}
  t := 0.0; {start at zero time}
  xold := 0.0; {start at zero position}
  x := xold; {and at zero velocity}
  d := 0.0; {assume no disturbance initially}
  v := 0.0;
  d1 := 1.5 + 6.0*random;
  d2 := d1 + 3.0;
  maxu := 0.0;
  legends;
  str(mslope:5:3,numstr);
  setcolor(white);
  outtextxy(0,60,'Slope = '+numstr);
  outtextxy(150,100,'MCT MODEL');

repeat {the control loop starts here}

  r := make_reference (t+dt,mslope); {define reference}

  {COMPUTE OUTPUT}
  u := K * (r - 2.0 * x + xold) - d;

  {LIMIT OUTPUT TO +/- 100 N-M}
  if u < -100.0 then u := -100.0 else {limit output, if desired}
    if u > +100.0 then u := +100.0;

  {GENERATE PREDICTED X}
  xpre := 2.0 * x - xold + u / K + d / K;

  {ENVIRONMENTAL EQUATIONS}
  do_observation;
  plotit(110);
  t := t + dt;

  {ESTIMATE DISTURBANCE}
  d := d + K * (x - xpre);

until t >= 9.0; {at this point the loop ends}

ch := readkey;
if ch in ['=','+'] then mslope := mslope*1.5;
if ch in ['_','-'] then mslope := mslope/1.5;

until ch in ['q','Q','p','P'];
end;

{initialization}
begin
  clrscr; {clear screen}
  initgraphics;
  maxy := getmaxy;
  maxx := getmaxx;
  ycenter := (maxy+1) div 2;
  J := 10.0; {or whatever value...}
  dt := 0.01; {or whatever value...}
  K := (J / dt) / dt; {auxiliary constant}
  mslope := 0.2;
  pslope := 0.2;
  setcolor(white);
  outtextxy(0,maxy - 15,'q to quit, space to repeat, +/- to change slope');
  ch := 'm';
  repeat
   if ch in ['m','M'] then MCTmodel;
   if ch in ['p','P'] then PCTmodel;
  until ch in ['q','Q'];
  closegraph;
end.

···

a := (u+trued)/J;

[HANS_BLOM2]

[Hans Blom, 970306]

(Bill Powers (970305.0501 MST))

Bill, this nonsense about "exact" has got to stop.

The physically correct environmental model is a simple application
of the laws of motion. If a constant torque u is applied to the
theodolite for a time dt, the angle x will vary as

x := x0 + v*dt + 0.5*a*dt*dt

This is a _program step_, not a mathematical relationship. In this
program step, times or time indices have been left out. What is the
mathematical meaning of your line? Assuming that x0 belongs to time
t, does x belong to time t as well? No, obviously not; in that case,
the line would just be x := x0. So x belongs to time t+dt. Let's plug
this in:

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

Now we have to answer the questions which times go with v and a. Do
we need to write v(t) or v(t+dt), and a(t) or a(t+dt)? Or v(t+0.5*dt)
maybe? In the program (and in the difference equation) continuous
time is chunked into intervals of size dt. So v(t+0.5*dt) is not
available, regrettably, and that does not depend on the choice of dt;
it is either v(t) or v(t+dt). The same with a, of course.

Does it matter? Yes: v and a will generally vary continuously
throughout the interval dt, even though that variation may be small
(with a small dt). So not any value of q, 0<x<1, in v(t+q*dt) will
do; different values of q will give different outcomes. If the output
is to be the same regardless of the value of q, this can only be done
by postulating that x, v and a are constant throughout an interval
dt. Our initial conditions and the program take care of this being
true for a, but not for v and x. So the only proper interpretation of
your program line is

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

which is identical to my prediction equation under the condition that
v(t) is known. And exactly that is the problem!

What all this means is: in the program, and in the difference
equation(s), the values of x and v are only known at a number of
discrete time instants T, T+dt, T+2dt, ... At all other times the
values are unknown and cannot be talked about (within this
formalism). If you do, you just generate confusion. Or you assume
something, e.g. that the values do not change during the interval.
The latter is incorrect, even though it may be a good approximation
in practice, particularly if dt is small. But it is not _exact_ and
it will introduce what you have called "integration errors" in your
approach.

The same discussion with the next equation:

The velocity v will increase as

v := v + a*dt.

Its proper interpretation is

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

The angular acceleration a is just (u + d)/J: torque divided by
moment of inertia. Since (u + d) is constant during the interval dt,
this is the exact way to calculate the value of x at the start of
the next iteration.

In formula, this is

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

This is not a differential but an algebraic equation: all terms refer
to the same time.

Why is all this important? It becomes important when x and v are not
just hypothetical things (you treat v that way!) but variables that
we need to _perceive_. How do we perceive v(t)? Our sensors only
provide us with values for x. Many estimates are possible for v(t) at
time t (but remember that perceptions _beyond_ time t are not
allowed!), such as

v(t) = [x(t) - x(t-dt)]/dt
v(t) = [x(t) - x(t-2dt)]/2dt
v(t) = [x(t) - x(t-dt) + x(t-2dt) - x(t-3dt)]/2dt

or any other "perceptual input function" that you might want to
define. But in all cases will such a constructed "perception" of v(t)
run behind its real physical equivalent; in the first expression, for
example, by 0.5dt, and in the second one by 1.0dt. This running-
behind approximation is inherent in an _observer_; it does not
describe the _physical_ model. [By the way: it introduces a
Heisenberg-like uncertainty relation if the observation is noisy or
obtained with finite precision or resolution.] Yet you use it for
that purpose. That is not _exact_, as you pretend it is.

This being the case, both the MCT and the PCT models should use the
exact equation of motion, and not the approximation in the original
MCT model.

I maintain that the expression that I used in my MCT program is an
exact translation of the original second order differential equation
J*d2/dx2=u+d and independent of dt (with the only assumption that the
system's acceleration u+d changes at the boundaries of our time
intervals and not within them), whereas your approximation is not.

How to test that? In the PCT model, the open-loop "PCT-physical" x
will run slightly behind the "true physical" x. The latter,
regrettably, is not available -- unless one uses a tool such as
Matlab which can accurately simulate differential equations. We know,
however, that the running-behind depends on the choice of dt. So at
different values of dt, particularly rather large ones (e.g. dt=0.1
or dt=1.0), the behavior of your (open loop!) environment equation
(say with a constant a=1) will be dependent on the value of dt.
Whereas this is not so for the equation I used. I therefore maintain
that my equation is exact (given the assumption stated above) and not
contaminated by dt-truncation effects, whereas yours is. Please
check...

An adjustment has to be made to the "predictor" equation to correct
for the problems that come from using the exact environmental
equations in the MCT model. I will leave it to Hans to make the
adjustment.

An adjustment has to be made to the environment equation to correct
for its inexactness in the PCT model. I will leave it to Bill to make
the adjustment. Anyway, we started this whole discussion with
granting the MCT a perfect model, i.e. a perfect match between the
true environment function and its internal representation of the
environment. If you change the environment function, the model isn't
perfect anymore.

Your exercise does have a positive side, however: it shows how this
MCT controller operates if the model is _not_ perfect.

By the way, your new program still has the

  j := 1.0

which means that the MCT controller must work with a ten times
heavier load than the PCT controller. That's not fair!

I'm sorry to see that we remain stuck in side issues for so long.

Greetings,

Hans

[HANS_BLOM2]

[From Bill Powers (970305.0501 MST)]

Attached is theo5mct.pas, which is a rearrangement of theo4mct.pas.

Nice work, Bill. I changed the code such that J is defined as 1 for
both controllers. As I expected, their behavior is almost identical.
The MCT controller is more precise, reaches the reference level
slightly sooner on very fast ramps, it seems, but at the cost of more
overshoot. Your controller is, indeed, well tuned! What's your
secret? Trial and error?

It is only when different Js need to be controlled -- which needs
changing the constant that defines J and recompilation of the program
now -- that the MCT controller convincingly shows its superiority.
But that is no surprise, of course. It knows J, after all...

In "do_observation," it will be seen that there are _two_
environmental equations, one used for the MCT model and the other
for the PCT model.

Yes, and they are slightly different, as you note. But I discussed
that already in my earlier post.

Appended is the new program; I will delay posting it with the
runnable version to my FTP page until Hans has the MCT model
running with the correct environmental model.

The MCT model _has_ the correct environmental model!

If I cannot convince you of that, I'm willing to use your approach as
the "real" thing -- even though it's not -- and design my controller
with a model that matches it, again perfectly. But I'd much rather
not do that, of course. It would achieve nothing except that these
equations can be "inverted" as well. Moreover, the pile of work on my
desk is unusually high at the moment...

Greetings,

Hans

[Bill_Powers1]

[From Bill Powers (970306.0632 MST)]

Hans Blom, 970306--

Bill, this nonsense about "exact" has got to stop.

Not bloody likely! :^)>

The physically correct environmental model is a simple application
of the laws of motion. If a constant torque u is applied to the
theodolite for a time dt, the angle x will vary as

x := x0 + v*dt + 0.5*a*dt*dt

This is a _program step_, not a mathematical relationship. In this
program step, times or time indices have been left out. What is the
mathematical meaning of your line? Assuming that x0 belongs to time
t, does x belong to time t as well? No, obviously not; in that case,
the line would just be x := x0. So x belongs to time t+dt.

Yes, and x0 is x(t), where t is a _continuous variable_. Consider the output
of a discrete controller applied as a torque to our theodolite:
                                              .
                                             .
                                           . ----------------------
                                         . |
      ^ -----------------.-----
      > > .
     u> > . <----x
-------------------- . . . .
t1 t2 t3 t4
                      <-------delta-t ------->

This discrete output is applied to a continuous device, which reacts to
torques as

x = x0 + a*dt + 0.5*a*dt^2

where dt is interpreted as an infinitesimal time interval carried to the
limit of zero. The behavior of x is approximated above by the dots.

Assume that the initial velocity and torque shown above are zero, at t1.
The angle will remain at zero until t2. Then we can say

x(t3) = x(t2) + v(t2)*dt + 0.5*a(t2)*dt^2, where

v(t2) = v(t1) + a(t1)*dt, and

a(t) = u(t)/J.

This is the EXACT response of the physical theodolite to the stepped torque
input, under the usual assumptions (no friction, noise, etc.).

When, inside the MCT system, you simulate this response by using discrete
equations, you replace the parabolic segments of the actual response by
straight-line approximations, and the infinitesimal dt by a finite delta-t.
As we all learned back in introductory calculus, this introduces an error.
The error gets smaller as delta-t approaches zero, but it never becomes zero
until delta-t = 0.

What all this means is: in the program, and in the difference
equation(s), the values of x and v are only known at a number of
discrete time instants T, T+dt, T+2dt, ... At all other times the
values are unknown and cannot be talked about (within this
formalism). If you do, you just generate confusion.

It's NOT talking about them that generates confusion. The values of the
angle and velocity are known _to the discrete controller_ only at those
times, but to an external observer working with physical instruments, the
continuous motion of the theodolite is the basic observation. We can measure
to milliradians with a resolution of microseconds if need be; the motion
will appear continuous on any scale of observation larger than that of
elementary particles. The physical theodolite does not behave differently
simply because it is observed only at finite intervals, or subject to
stepped forces instead of continuously-variables ones.

Or you assume
something, e.g. that the values do not change during the interval.
The latter is incorrect, even though it may be a good approximation
in practice, particularly if dt is small. But it is not _exact_ and
it will introduce what you have called "integration errors" in your
approach.

That's an _additional_ source of error, a random error. What we are talking
about here is a _systematic_ error introduced by the discrete approximation
to a continuous process, even if the output is constant between steps.

I maintain that the expression that I used in my MCT program is an
exact translation of the original second order differential equation
J*d2/dx2=u+d and independent of dt (with the only assumption that the
system's acceleration u+d changes at the boundaries of our time
intervals and not within them), whereas your approximation is not.

If you actually solve that differential equation using the calculus, you
will find that the solution is exactly in the form I state. As I have shown
above, because of the stepped nature of the discrete controller's output the
response of the physical theodolite will be a series of parabolic segments,
not straight lines. I appeal to our mathematical readers: is this not true?

The test of the equivalence of your expression to mine is to substitute mine
for yours in your model. This is easily done in theo5mct, by inserting a
"not" before "ismct" in the "do_observation" procedure (formerly a
function). You will find that your model immediately begins to show strong
oscillations. The only reason it doesn't oscillate with your model of the
theodolite is that you have assumed your approximation to apply to the
theodolite itself, so the theodolite actually moves in straight-line
segments during each dt. This FALSE ASSUMPTION is what makes your model work.

Your analysis depends on there being a mathematically EXACT correspondence
of the controller's simulation of the theodolite to the actual response of
the theodolite to the stepped output. If there is even a tiny difference,
your model will attempt to correct it, and in doing so will introduce an
even larger difference, the end-result being that that output alternates
back and forth between maximum positive and negative values. Try it and see.

It is possible that some small adjustment in your model can eliminate this
problem. I don't know; it's up to you to find it. I'm not going to make a
prediction one way or the other.

I highly recommend doing the comparisons you mention, using a precision
integrator like that in Matlab. I compared your model with mine using a
constant applied torque, and indeed found differences in the output, ranging
up to 4% after 1 second, with dt = 0.01; the difference gets smaller (as
expected) as dt decreases. I believe that the Matlab simulation will support
my exact solution of the differential equation.

How to test that? In the PCT model, the open-loop "PCT-physical" x
will run slightly behind the "true physical" x.

We aren't talking about the PCT model, but the physical model of the
theodolite. This is just a simple matter of elementary physics, not of
differences between controller designs. Both of our models apply stepped
output torques to the theodolite, and both sample the angle -- at the same
intervals. All I am doing is _correctly_ computing the response of the
theodolite to the torque steps, independently of where the torques are
coming from. Your model computes it _incorrectly_ -- not by much, but
apparently by a critical amount when it concerns how your model functions.

By the way, your new program still has the

j := 1.0

which means that the MCT controller must work with a ten times
heavier load than the PCT controller. That's not fair!

I did correct it, and said that after the correction the two models behaved
essentially identically. I haven't posted that corrected code, but it's
appended below.

The peak torques are now shown at the left of the screen. It's interesting
that at the same slope (now initialized to 0.4), varying dt has opposite
effects on the peak output torques:

  dt MCT max u PCT max u
  0.01 sec 150 n-m 72 n-m
  0.002 200 62
  0.0005 2900 61
  0.0001 78099 60

As dt goes toward zero, the peak MCT torque goes toward infinity; the peak
PCT torque tends toward a minimum. The MCT model also oscillates severely at
the start and the end of the ramp when dt gets very small.

I'm sorry to see that we remain stuck in side issues for so long.

Perhaps you'll come to see them as more important. I hope the appended code
doesn't contain any more of my mistakes.

Best,

Bill P.

···

=======================================================================
program theo5mct;
{
Mistaken initialization corrected 970306 WTP. Substitute for old
version.

This program compares the performance of the MCT and PCT models. Both
models are controlling the position of a theodolite, under the
following conditions:

change of position: 1 radian.
time resolution: 0.01 sec (dt).
Maximum available output torque: 100 newton-meters
Moment of inertia of theodolite: 1 newton-meter^2
Rate of change of reference signal: variable from keyboard
disturbance: 100 n-m, 3 sec duration, random starting time.

The MCT model runs first, in the lower half of the screen. The slope
of the change in reference signal can be adjusted up and down (by
factors of 1.5) by pressing the + and - keys, shifted or unshifted. A
new run is done every time any key is pressed -- except 'p' or 'q'.

Pressing the 'q' key exits the program at any time. Pressing the 'p'
key switches to running the PCT model; after the first time, pressing
'p' or 'm' switches control to the PCT or MCT model, for further
adjustments of the slopes.

The beginning time of the disturbance varies randomly each time a new
run occurs. By pressing the space bar, new runs can be done with the
same value of the slope of the reference signal's rise but a different
start of the disturbance. Thus you can see what happens when there
is or is not a disturbance occurring during the change in reference
signal.

MCT program by Hans Blom; PCT program and presentation by W. T. Powers
4 March 1997
}

uses
  dos, crt, graph, grutils;

var
  J, K, dt, r, x, u, a, v, xold, xpre, xsav: real;
  d, d1, d2, trued, t, pslope,mslope, umax: real;
  gv, gx, rv, rx: real;
  xplot: integer;
  maxx, maxy,ycenter: integer;
  ch: char;
  numstr: string;
  fo: text;
  ismct: boolean;

function make_reference (t: real; var slope: real): real;
  {this function defines the setpoint at time t}
var ref: real;
begin
if (t >= 3.0) then
begin
  ref := slope*(t - 3.0);
  if ref > 1.0 then ref := 1.0;
  make_reference := ref;
end
  else make_reference := 0.0
end;

function true_disturb (t: real): real;
  {this function defines the true disturbance}
begin
  if (t < 6.0) or (t > 8.0) then
    true_disturb := 0.0
  else
    true_disturb := 50.0;
end;

procedure do_observation;
{this function generates the x that the controller will observe}
begin
  trued := true_disturb(t);
  xsav := x; {save present x}
  if ismct then
   begin
    x := 2.0 * x - xold + u/K + trued/K;
   end
  else
   begin
    a := (u+trued)/J;
    x := x + v*dt + 0.5*a*dt*dt;
    v := v + a*dt;
   end;
  xold := xsav; {xold := previous x}
end;

procedure plotit(baseline: integer);
begin
  xplot := round(50.0*t)+ 180;
  putpixel(xplot,baseline,white);
  putpixel(xplot,baseline - round(100.0*r),yellow);
  putpixel(xplot,baseline - round(100.0*x),white);
  putpixel(xplot,baseline - round(trued),lightred);
  putpixel(xplot,baseline - round(u),lightcyan);
end;

procedure legends;
begin
  clearviewport;
  setcolor(white);
  outtextxy(0,0,'Pointing angle');
  setcolor(yellow);
  outtextxy(0,15,'Ref level');
  setcolor(lightcyan);
  outtextxy(0,30,'Output torque');
  setcolor(lightred);
  outtextxy(0,45,'Disturbance');

end;

Procedure PCTmodel;
begin
ismct := false;
setviewport(0,0,maxx,ycenter,ClipOn);
repeat {REPEAT WHOLE RUN}
  t := 0.0; {start at zero time}
  d := 0.0; {assume no disturbance initially}
  v := 0.0;
  x := 0.0;
  gv := 100.0;
  gx := 50.0;

  legends;
  setcolor(white);
  line(0,ycenter-8,maxx,ycenter-8);
  str(pslope:5:3,numstr);
  setcolor(white);
  outtextxy(0,60,'Slope = '+numstr);

  outtextxy(150,100,'PCT MODEL');
  umax := 0.0;
repeat {the control loop starts here}

  r := make_reference (t,pslope); {define reference}

{CONTROL SYSTEM EQUATIONS}
  rv := gx* (r - x); {velocity ref level = output of position control}
  u := gv*(rv - v); {output force = output of velocity control}

  if abs(u) > umax then umax := abs(u);

  if u < -100.0 then u := -100.0 else {limit output, if desired}
    if u > +100.0 then u := +100.0;

  {ENVIRONMENTAL EQUATIONS}
  do_observation;
  plotit(120);
  t := t + dt;

until t >= 9.0; {at this point the loop ends}
str(umax:1:0,numstr);
outtextxy(0,75,'Max torque = ' + numstr + ' n-m');

ch := readkey;
if ch in ['=','+'] then pslope := pslope*1.5;
if ch in ['_','-'] then pslope := pslope/1.5;
until ch in ['q','Q','m','M'];
end;

Procedure MCTmodel;
begin
ismct := true;
setviewport(0,ycenter+1,maxx,maxy - 10,ClipOff);
repeat {REPEAT WHOLE RUN}
  t := 0.0; {start at zero time}
  xold := 0.0; {start at zero position}
  x := xold; {and at zero velocity}
  d := 0.0; {assume no disturbance initially}
  v := 0.0;
  legends;
  str(mslope:5:3,numstr);
  setcolor(white);
  outtextxy(0,60,'Slope = '+numstr);
  outtextxy(150,100,'MCT MODEL');
  umax := 0.0;
repeat {the control loop starts here}

  r := make_reference (t,mslope); {define reference}

  {COMPUTE OUTPUT}
  u := K * (r - 2.0 * x + xold) - d;
  if abs(u) > umax then umax := abs(u);

  {LIMIT OUTPUT TO +/- 100 N-M}
  if u < -100.0 then u := -100.0 else {limit output, if desired}
    if u > +100.0 then u := +100.0;

  {GENERATE PREDICTED X}
  xpre := (2.0 * x - xold + u / K + d / K);

  {ENVIRONMENTAL EQUATIONS}
  do_observation;
  plotit(125);
  t := t + dt;

  {ESTIMATE DISTURBANCE}
  d := d + K * (x - xpre);

until t >= 9.0; {at this point the loop ends}

str(umax:1:0,numstr);
outtextxy(0,75,'Max torque = ' + numstr + ' n-m');

ch := readkey;
if ch in ['=','+'] then mslope := mslope*1.5;
if ch in ['_','-'] then mslope := mslope/1.5;

until ch in ['q','Q','p','P'];
end;

{initialization}
begin
  clrscr; {clear screen}
  initgraphics;
  maxy := getmaxy;
  maxx := getmaxx;
  ycenter := (maxy+1) div 2;
  J := 1.0; {or whatever value...}
  dt := 0.01; {or whatever value...}
  K := (J / dt) / dt; {auxiliary constant}
  mslope := 0.4;
  pslope := 0.4;
  setcolor(white);
  outtextxy(0,maxy - 10,'q to quit, space to repeat, +/- to change slope');
  ch := 'm';
  repeat

  if ch in ['m','M'] then MCTmodel;
  if ch in ['p','P'] then PCTmodel;
  until ch in ['q','Q'];
  closegraph;
end.

[HANS_BLOM2]

[Hans Blom, 970310]

(Bill Powers (970306.0632 MST))

The test of the equivalence of your expression to mine is to
substitute mine for yours in your model. This is easily done in
theo5mct, by inserting a "not" before "ismct" in the
"do_observation" procedure (formerly a function). You will find
that your model immediately begins to show strong oscillations.

Is there a better demonstration that the two expressions are _not_
equivalent? They may look that way, and your analysis may say that
they are the same, but when, on exchange of one by the other,
different behavior results, they are obviously not the same.

Now the question: what is the difference?

Greetings,

Hans