# Revised theodolite model

[From Bill Powers (970226.0945 MST)]

Can't leave well enough alone. At the end is the source code for a new
version of the PCT theodolite program; I'll also substitute the .exe file on
my ftp page, at www.frontier.net/users/powers_w. Just replace the old copies
with the new ones.

This version allows you to change dt from the keyboard. The upper limit is
set so you will see some computational oscillations when dt gets too large
for the gains of the control loops. There is no lower limit, so if you've
got an hour or two you can set dt to 1e-10 if you want. The + and - keys
double and halve the value of dt. Striking any key cuts the run short, so
you can raise and lower dt quickly. Also, the legends now include the
instructions (!) -- q to quit, space for another run, any key during the run
cuts it short.

Best,

Bill P.

···

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

{corrected for long integer in "for" loop to allow shorter dt}
{instructions added to screen}
{allows changing dt from keyboard}

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: longint;
k,xscale: integer;
ch: char;
numstr: string;
begin
initgraphics;
dt := 0.001; {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}
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}
xscale := round(9.0/dt/450.0); {x scaling for plots}
clearviewport;
setcolor(white);
outtextxy(0,0,'Pointing angle');
setcolor(lightcyan);
outtextxy(0,15,'Output torque');
setcolor(lightred);
outtextxy(0,30,'Disturbance');
setcolor(white);
str(dt:8:6,numstr);
outtextxy(0,45,'dt = ' + numstr + ' sec');
outtextxy(0,240 - 85,'1 radian');
outtextxy(0,240 - 104,'---------->');
setcolor(lightcyan);
outtextxy(0,240 - 124,'400 newton-meters');
setcolor(white);
outtextxy(0,420,'q to quit, space to repeat, any key to cut short');
outtextxy(0,435,'+,- to double or halve dt');

{ITERATION LOOP FOR RUN}
for i := 1 to round(9.0/dt) do
begin
if keypressed then i := round(9.0/dt);
{CHECK FOR CHANGES IN REF SIGNAL, DISTURANCE}
if (t > r1) and (t < r2)
then rx := 1 else rx := 0.0; {POSITION 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 +/- 500 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}
if ch in ['+','='] then dt := 2.0*dt;
if ch in ['-','_'] then dt := 0.5*dt;
if dt > 0.004 then dt := 0.004;

until ch = 'q'; {repeat on any key but 'q'}

closegraph;
end.

[Hans Blom, 970227c]

Bill, there's something that I don't understand in this program (or
the previous one, for that matter): where is the J? In other words:
how does the program mimick the response of the theodolite when u is
applied? What is the next x going to be in terms of J? In other words
again: what would need to be changed in the program if J happens to
be twice larger or smaller than you now assume?

Greetings,

Hans

[From Bill Powers (970227.0931 MST)]

Hans Blom, 970227c--

Bill, there's something that I don't understand in this program (or
the previous one, for that matter): where is the J? In other words:
how does the program mimick the response of the theodolite when u is
applied? What is the next x going to be in terms of J? In other words
again: what would need to be changed in the program if J happens to
be twice larger or smaller than you now assume?

Shoot, you caught it before I could send my correction off. See previous
post, and P.S. at the end. I'll send a corrected version soon.

Best,

Bill P.