PCT model of theodolite controller

[From Bill Powers (970225.0930 MST)]

Hans Blom, Rupert Young (970225) and any other interested programmers--

At the end is the Pascal source code for my PCT model of controlling
theodolite angle. I have added one new condition and changed one, after
seeing how the simulated theodolite behaves under applied forces. The new
condition is that the output torque is limited to +/- 750 newton-meters; the
changed condition is that the disturbance magnitude is 100 newton-meters
when it is present, instead of 10.

The reason for the torque limit is that when we present the system with a
step-change in the reference signal, the torque goes instantly to a very
large magnitude (if the system is set for fast control), like many thousands
of newton-meters. That would be OK, I suppose, except that scaling the plot
to show the peak output makes the response to the disturbance invisible.
I've compromised by limiting the output to 750 newton-meters (about 450
foot-pounds, which is like attaching a Chevy truck engine to turn the
theodolite) and increasing the disturbance to 100 newton-meters, to make the
response to it easily visible. I don't think these changes put any undue
limits on a realistic design.

I have not put any phase advance into the PCT model. It doesn't seem
necessary. Any more torque would probably bend the tube of a real theodolite.

It's easy to miss the initial response of the system. The output, in cyan or
light blue, jumps to the maximum for a few dt's when the reference signal
jumps to 1 radian, making a dot, and then falls to the maximum negative or
near it, approaching zero from below as it decelerates the theodolite to its
final position. When the reference signal drops back to zero, the reverse
sequence takes place. Note that when the theodolite is stationary (white
line horizontal), the output torque is zero except when the disturbance is
present.

The duration of the run is 9 seconds.

Feel free to incorporate any or all of this code into your programs. It
requires the Borland .bgi files to be in the same directory where the
program runs. If you don't have the "grutils" unit already, I will post it
again.

Best,

Bill P.

···

========================================================================
program theodPCT;

uses dos, crt, graph, grutils;

var j,u,x,rx,v,rv,d,t,dt,dv,dx,gv,gx,r1,r2,d1,d2: real;
    i,k,xscale: integer;
    ch: char;
begin
initgraphics;
dt := 0.002; {time resolution}
r1 := 3.0; {start of ref signal change}
r2 := 6.0; {end of ref signal change}
gv := 500.0; {gain of position control loop}
gx := 40.0; {gain of velocity control loop}
xscale := round(9.0/dt/450.0); {x scaling for plots}
randomize;
repeat
rv := 0.0; {initialize velocity ref level & other variables}
x := 0.0;
v := 0.0;
t := 0.0;
d1 := 1.5 + 6.0*random; {start of disturbance}
d2 := d1 + 2.0; {end of disturbance}
clearviewport;
setcolor(white);
outtextxy(0,0,'Pointing angle');
setcolor(lightcyan);
outtextxy(0,15,'Output torque');
setcolor(lightred);
outtextxy(0,30,'Disturbance');
setcolor(white);
outtextxy(0,240 - 85,'1 radian');
outtextxy(0,240 - 104,'---------->');
outtextxy(0,240 - 120,'400 newton-meters');

{ITERATION LOOP FOR RUN}
for i := 1 to round(9.0/dt) do
begin

{CHECK FOR CHANGES IN REF SIGNAL, DISTURANCE}
if (t > r1) and (t < r2)
  then rx := 1 else rx := 0.0; {REF SIGNAL, RADIANS}

if (t > d1) and (t < d2) then
  d := 100 else d := 0; {DISTURBANCE, NEWTON-METERS}

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

{LIMITS ON OUTPUT}
if u > 750.0 then u := 750.0 else {limit output to +/- 750 n-m}
if u < -750.0 then u := -750.0;

{ENVIRONMENT EQUATIONS}
v := v + (u + d)*dt; {integrate torque to get angular velocity}
x := x + v*dt; {integrate angular velocity to get angle}

{PLOT VARIABLES}
putpixel(i div xscale, 240,white);
putpixel(i div xscale, 240 - round(x*100),white);
putpixel(i div xscale, 240 - round(d/4),lightred);
putpixel(i div xscale, 240 - round(u/4),lightcyan);

{ADVANCE TIME BY DT}
t := t + dt;

end; {END OF RUN LOOP}

ch := readkey; {pause to view}
until ch = 'q'; {repeat on any key but 'q'}

closegraph;
end.