[From Bill Powers (980225.0806 MST)]

Richard Kennaway (980225.2045 JST)--

Ah, you're safely in Japan. I know your thoughts are going to be mainly on

other subjects for the next six weeks, but perhaps when you need a break

you can think a little about PCT.

I put together a simulation of a pendulum, to serve as a test-bed for an

upside-down pendulum balancing control system (seems that everybody has to

try this once). This brought back an approach to simulation that I tried

with integer arithmetic some years ago, but now can try in floating point

because of a faster computer. I'd like to know what you think of it, and

whether you think it can be generalized to simulating mechanical systems of

other kinds (like legs and arms and bodies).

The pendulum is mounted on a cart at the upper end, and hangs below it on a

rod with a bob on the end. The rod is unbendable but slightly elastic along

its length. The cart has a mass of 0.1 kg, and the bob weighs 1.0 kg; the

rod is massless for now.

It's the elasticity of the rod that makes this different from other

simulations I've seen (you might like to think about how you'd do this, to

compare with the simulation below). The position of the bob is the position

of the mounting point on the cart plus the position of the bob relative to

the mounting point, in x and y (this is two-dimensional for now). The force

acting on the bob is aligned with the direction of the rod, and is equal to

(L - Lo)*ke newtons,

where Lo is the resting length, L is the actual length, and ke is the

constant of elasticity in newtons per meter.

This force is resolved into an x and a y component, and applied to the bob.

Only similar triangles are used -- no trigonometry at all.

The x and y forces produce x and y accelerations of the bob, and the

opposite acceleration of the cart in x (the cart is held on rails in y).

Gravity acts downward on the bob, too. Integrating these forces twice gives

the velocity and position of the bob, and the x position of the cart.

As you will see if you can compile and run the program below, the result

looks very realistic; it might even be physically correct. The actual

dynamical simulation is in the procedure "dynamics."

Note that the rod is very stiff: its coefficient of elasticity is a million

newtons per meter of stretch, or a thousand newtons per millimeter, or one

newton per micron. This approximates the proverbial "stiff rod," but allows

calculating force as a function of rod elongation. It is more physically

correct than the idealized rod, and seems to allow a nice shortcut in the

simulation. Right now I can't think of how I'd do this without the elastic

rod -- I suppose I've have to use centrifugal forces and so on. That's why

I'd like to know how you would do it on the basis of what's in those

classical mechanics books. I think this way is probably a lot simpler.

When the program comes up, the bob is held deflected to one side. Pressing

a key releases the bob and starts the run. Pressing a key again exits the

program. The red trace shows the cart position as it is dragged back and

forth by the pendulum swinging from it. There is a coefficient of friction

that exerts a braking force on the cart proportional to its velocity, so

the oscillations gradually damp out.

There's a very interesting phenomenon that occurs. The bob swings back and

forth, and as the energy of the system is dissipated by the friction in the

cart, the oscillations of the bob die out until it's hanging straight down

from the cart. At the same time, the cart oscillates back and forth,

acquiring a mean velocity which also eventually dies out to zero. What's

interesting is that the cart returns exactly to its initial position before

the bob is released! My intuitive guess would have been that the final

position of the cart would be the mid-point of its oscillations, but it is

not. This doesn't depend on the relative masses of the cart or the bob, or

the coeffient of friction, or the amplitude of the initial displacement of

the bob.

When I get to the control system, it will push on the cart, also through a

stiff spring. I wish we could create force feedback to the mouse, because

then you'd feel exactly how this dynamic system would feel to a hand

pushing on it. But it should still be possible to balance the pendulum

upside down by operating the mouse, and then compare that performance

(after learning is complete) with a simulated control system doing the same

thing.

This is another answer to the skeptics who think that PCT applies only to

tracking and linear systems. Of course they won't remember it.

Best,

Bill P.

## ···

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

program pendcart;

uses crt,dos,graph,grutils;

var bob: record x, y, dx, dy, fx, fy, vx, vy, mass: real end;

cart: record x, y, dx, dy, fx, fy, vx, vy, mass: real end;

dt: real;

l,l0,f,g,ke: real;

iters: longint;

ch: char;

procedure initpendulum;

begin

dt := 2e-4; { sec }

ke := 1e6; { newton/meter}

l0 := 1.0; { meter}

g := 9.8; { meter/sec^2}

iters := 0;

with bob do

begin

vx := 0.0; {m/sec}

vy := 0.0;

x := -l0*sin(0.8); {l0 = L-nought, not 10.0}

y := -l0*cos(0.8); {meter}

mass := 1.0 {kg}

end;

with cart do

begin

vx := 0.0;

vy := 0.0;

x := 0.0;

y := 0.0;

mass := 0.1;

end;

end;

procedure showpendulum;

const oldx: integer = 0;

oldy: integer = 0;

oldcx: integer = 0;

oldcy: integer = 0;

begin

if iters > 0 then

begin

setcolor(black);

rectangle(hcenter - 30 + oldcx,

vcenter,hcenter + 30 + oldcx,vcenter- 15);

circle(hcenter + oldx, vcenter - oldy,10);

circle(hcenter + oldcx - 20, vcenter - oldcy,7);

circle(hcenter + oldcx + 20, vcenter - oldcy,7);

line(hcenter + oldcx,vcenter - oldcy,

hcenter + oldx,vcenter - oldy);

end;

oldx := round(200*bob.x);

oldy := round(200*bob.y);

oldcx := round(200*cart.x);

oldcy := round(200*cart.y);

setcolor(white);

rectangle(hcenter - 30 + oldcx,

vcenter,hcenter + 30 + oldcx,vcenter- 15);

circle(hcenter + oldx, vcenter - oldy,10);

circle(hcenter + oldcx - 20, vcenter - oldcy,7);

circle(hcenter + oldcx + 20, vcenter - oldcy,7);

line(hcenter + oldcx - 20,vcenter + 7,hcenter + oldcx + 20,vcenter+ 7);

line(hcenter + oldcx,vcenter - oldcy,

hcenter + oldx,vcenter - oldy);

line(0, vcenter + 7,hsize,vcenter + 7);

end;

procedure dynamics;

begin

with bob do

begin

dx := bob.x - cart.x;

dy := bob.y - cart.y;

l := sqrt(sqr(dx) + sqr(dy));

f := (l0 - l)*ke; {l0 = L-nought, not 10.0 -- sorry}

fx := f*dx/l;

fy := f*dy/l - mass*g;

vx := vx + fx/mass* dt;

x := x + (vx + 0.5*fx*dt/mass)*dt;

vy := vy + fy* dt;

y := y + (vy + 0.5*fy*dt/mass)*dt;

end;

with cart do

begin

fx := -bob.fx - 0.1*vx;

vx := vx + fx/mass*dt;

x := x + (vx + 0.5*fx/mass*dt)*dt;

end;

end;

begin

initgraphics;

initpendulum;

l := l0; {l0 = L-nought, not 10.0 -- sorry}

outtextxy(0,0,'PRESS ANY KEY TO START OR QUIT');

setcolor(lightred);

outtextxy(0,50,'CART POSITION');

setcolor(white);

while not keypressed do

begin

dynamics;

if (iters mod 50) = 0 then

begin

showpendulum;

putpixel(iters div 500, 70 - round(100*cart.x),lightred);

putpixel(iters div 500, 70,white);

end;

if iters = 0 then

begin

line(hcenter,vcenter- 20,hcenter,vcenter - 35);

ch := readkey;

end;

inc(iters);

end;

closegraph;

end.