CONTROLLING THE INVERTED PENDULUM:
An Example of a Digital Feedback Control System.

A free swinging pendulum is inverted so the hinge is at the bottom. The challenge is to move the bottom of the pendulum back and forth to keep the pendulum from falling over.

The hinge at the bottom of the pendulum is attached to a cart which moves back and forth on a track. The cart is connected to wires which are connected to an electric motor. The voltage applied to the electric motor is controlled by a computer. The observable properties are the current pendulum's angle from vertical and the current cart position.

Below is a diagram of the setup:

cart

This paper derives all the characteristic equations of motion for an inverted pendulum, all the transfer functions, Bode plots, state space representations, and concludes with an example of controlling a real setup.

After mastering this 2-dimensional problem you can expand this to a 3-dimensional problem. If you can master the task of controlling an inverted pendulum free to fall in any direction, then you can launch a rocket into outer space, controlling the thrust angle at the bottom so it doesn't fall over.



THE INVERTED PENDULUM AND CART SYSTEM


EQUATIONS OF MOTION: The Inverted Pendulum


p0_1


cg  =  Center of gravity
p  =  Pivot point
 
xsubcg  =  X coorordinate of center of gravity
 
xdot1cg  =  d1xcgdt First derivative of  xsubcg
 
xdot2cg  =  d2xcgdt Second derivative of  xsubcg
 
ysubcg  =  Y coorordinate of center of gravity
 
ydot1cg  =  d1ycgdt First derivative of  ysubcg
 
ydot2cg  =  d2ycgdt Second derivative of  ysubcg
 
xsubp  =  X coordinate of pivot point.
Note positive xsubp is toward right.
 
xdot1p  =  d1xpdt First derivative of  xsubp
 
xdot2p  =  d2xpdt Second derivative of  xsubp
 
ysubp  =  Y coordinate of pivot point.
 
ydot1p  =  d1ypdt First derivative of  ysubp
 
ydot2p  =  d2ypdt Second derivative of  ysubp
H  =  Force in Horizontal direction.
V  =  Force in Vertical direction.
l  =  1/2 length of pendulum
m  =  mass of pendulum
theta  =  Angle of pendulum measured from vertical.
Note positive theta is taken as clockwise.
bsubr  =  Viscous damping constant.
I  =  moment of inertia for pendulum     ( mip for uniform rod)
EQUATIONS OF MOTION: The Inverted Pendulum


Coordinates of point cg in terms of point p, angle theta, and length l

          p1a1

          p1a2

          p1a3


          p1a4

          p1a5

          p1a6



Sum of forces in X direction


          p1b1

          p1b2

          p1b3



Sum of forces in Y direction


          p1c1

          p1c2

          p1c3



Sum of moments about center of gravity


          p1d1

          p1d2



Combine Equations

          p2_1

          p2_2

          p2_3

          p2_4

          p2_5

Assume pendulum is uniform rod with moment of inertia  mip

Also assume small angle approximations sin( theta ) = theta , cos( theta ) = 1


          p2_6

          p2_7

Let   p2_8a,       p2_8b,       p2_8c


          p2_9

This is our characteristic equation of motion for the inverted pendulum.
TRANSFER FUNCTIONS

In control theory, functions called "transfer functions" are very often used to characterize the input-output relationships of linear time-invariant systems. The concept of transfer functions applies only to linear time-invariant systems, although it can be extended to certain nonlinear control systems.

Transfer Functions. The transfer function of a linear time-invariant system is defined to be the ratio of the Laplace transform of the output (response function) to the Laplace transform of the input (driving function), under the assumption that all initial conditions are zero.

The transfer function is an expression relating the output and input of a linear time-invariant system in terms of the system parameters and is a property of the system itself, independent of the driving function. The transfer function includes the units necessary to relate the input to the output; however, it does not provide any information concerning the physical structure of the system. (The transfer functions of many physically different systems can be identical.)

By using this concept, one can represent the system dynamics by algebraic equations in S.
          Modern Control Engineering by Katsuhiko Ogata, pg. 72. (University of Minnesota.) (©1970)



Derivation of Transfer Function for the Inverted Pendulum

Starting with our characteristic equation from the previous page:

          p2_9

we take the Laplace Transform of both sides:

          p3_2

and do some rearranging:

          p3_3

          p3_4


This is our transfer function for the inverted pendulum.
EQUATIONS OF MOTION: The Free Hanging Pendulum


p4_1

cg  =  Center of gravity
p  =  Pivot point
 
xsubcg  =  X coorordinate of center of gravity
 
xdot1cg  =  d1xcgdt   First derivative of  xsubcg
 
xdot2cg  =  d2xcgdt   Second derivative of  xsubcg
 
ysubcg  =  Y coorordinate of center of gravity
 
ydot1cg  =  d1ycgdt   First derivative of  ysubcg
 
ydot2cg  =  d2ycgdt   Second derivative of  ysubcg
 
xsubp  =  X coordinate of pivot point.
Note positive xsubp is toward right.
 
xdot1p  =  d1xpdt   First derivative of  xsubp
 
xdot2p  =  d2xpdt   Second derivative of  xsubp
 
ysubp  =  Y coordinate of pivot point.
 
ydot1p  =  d1ypdt   First derivative of  ysubp
 
ydot2p  =  d2ypdt   Second derivative of  ysubp
 
H  =  Force in Horizontal direction.
V  =  Force in Vertical direction.
l  =  1/2 length of pendulum
m  =  mass of pendulum
theta  =  Angle of pendulum measured from vertical.
Note positive theta is taken as clockwise.
bsubr  =  Viscous damping constant.
I  =  moment of inertia for pendulum     ( mip for uniform rod)



Coordinates of point cg in terms of point p, angle theta , and length l


          p5a1

          p5a2

          p5a3


          p5a4

          p5a5

          p5a6



Sum of forces in X direction


          p1b1

          p5b2

          p5b3



Sum of forces in Y direction


          p1c1

          p5c2

          p5c3



Sum of moments about center of gravity


          p1d1

          p5d2



Combine Equations


          p6_1

          p6_2

          p6_3

          p6_4

          p6_5

Assume pendulum is uniform rod with moment of inertia theta
Also assume small angle approximation sin( theta ) = theta


          p6_6

          p6_7

Let   p6_8a,       p6_8b


          p6_9

This is our characteristic equation of motion for the free hanging pendulum. Compare it with the characteristic equation for the inverted pendulum.


Solution of characteristic equation of motion for the free hanging pendulum


          p6_9

Assume a solution of the form   p7_2a
p7_2b
p7_2c


           p7_3

           p7_4

           p7_5

           Let  p7_6

           p7_7a            p7_7b


thetat  =  p7_8
 
 =  p7_9
 
 =  p7_10
 
 =  p7_11
 
 =  p7_12
 

Here p7_asub1, p7_asub2, and B are constants dependant upon initial conditions. In general the response  thetat   will be damped sinusoidal oscillations with a frequency omegd2p  and damping envelope of  dampenvl


Determination of Pendulum Parameters
zeta, omegan, and omegad

Remove the pendulum assembly from the cart and turn it upside down, allowing the pendulum to hang freely. Start the pendulum swinging and record the position  theta on a chart recorder. From our previous analysis we expect the response to be something like:

sindamp

The response will be damped sinusoidal oscillations with a frequency omegd2p and a damping envelope of dampenvl

From your chart recording estimate the parameters  zeta, omegan,and omegad

When I performed this experiment I found the decay appeared to be linear rather than exponential. It appeared the pendulum lost a fixed amount of energy every swing regardless of how wide the swing was.

viscous4

My guess is at the end of each swing when the pendulum momentarily stops the friction changes from kinetic (moving) friction to static (non-moving) friction, and it takes a certain amount of energy to overcome this static friction. zeta was estimated to be insignificantly small (0.02 as I recall, and I think we set zeta to zero to simplify the equations later on and ignored the static friction energy loss problem).



EQUATIONS OF MOTION: The Cart and Motor


cart
vsubi  =  input voltage. Positive input voltage moves cart toward right.
vsubo  =  output voltage
E  =  EMF across motor
rsuba  =  Motor armature resistance
isuba  =  Motor armature current
tsubd  =  Torque developed by motor
F  =  Force transmitted to wire by motor
phi  =  Motor shaft angle
r  =  effective radius of motor shaft. ( p9_delta )
H  =  Force in Horizontal direction.
xsubp  =  X coordinate of cart. Note positive  xsubp is toward right.
ksuba  =  Amplifier gain constant
ksubm  =  A constant of the motor
M  =  mass of cart
bsubm  =  viscous damping of motor
bsubc  =  viscous damping of cart
isubm  =  inertia of motor
J  =  combined inertia of motor and cart. ( p9_jeq )
B  =  combined damping of motor and cart. ( p9_beq )
ksubc  =  Viscous damping constant. ( p9_kceq )
THE CART AND MOTOR


Basic Equations of Motor and Amplifier


          p10a1

          p10a2

          p10a3

          p10a4



Combine Equations


          p10b1

          p10b2




Calcuate Torque


          p10c1

          p10c2

          p10c3


Assume m<<M

          p10c5

          p10c6

          p10c7

          p10c8




Substitute back in


          p11a1



          p11a2



          p11a3




Observe Following Relations


          p11b1

          p11b2

          p11b3




Substitute in


          p11c1



          p11c2



Let     p9_jeq
 
p11c4
 
p11c5


          p11c6

This is our characteristic equation of motion for the cart and motor system.

Determination of Parameters for the cart and motor


Now take the Laplace Transform of both sides of the characteristic equation on the previous page and rearrange to get the transfer function for the cart and motor.

The transfer function for the cart and motor is:

          p11d4

          or

          p13_1b



We shall experimentally determine the unknown parameters ksubc and B. The standard method for doing this is to input sinusoids of various frequencies into the system and make a Bode plot of the resulting output. We expect the Bode plot to have a breakpoint at the frequency corresponding to the unknown parameter B. The equations which predict this response are derived below.


           Let    p13_2


Here A is the amplitude of the sine wave we are inputting to the system and omega is the frequency in radians of the oscillations.

Taking the Laplace transform of both sides gives:

          p13_3

Applying this as the input to the transfer function gives:

           p13_4a p13_4b
 
p13_5


Expanding by partial fractions:

          p13_6

          p13_7

          p13_8

          p13_9

Thus:

          p13_10

Taking the inverse Laplace transform gives:


           p13_11a p13_11b
 
p13_12


The last term dies away as times increases, so we are left with a steady state solution of:

          p13_13

For those interested: Repeat the above derivation for an input of   p13_14   Explain the differences in the results.
The Bode Plot


When we perform the preceding experiment the output we measure is a voltage vsubx(t) which is related to the cart position X(t) by the equation:



          p15_1a

          or

          p15_1b

Thus our steady state output for an input

          acoswt

          is

          p15_2

We now wish to make a Bode plot of the magnitude (in decibels) of the output voltage corresponding to the cart position versus the frequency (in radians) of the input. The magnitude in decibels of vsubx(t) is

           p15_3a p15_3b
 
p15_4b
 
p15_5b
 
p15_6b




The Bode plot thus consists of three parts. First there is a constant gain term of p15_7a This is a horizontal line on the Bode plot. Second, there is a pole at the origin with a gain of p15_7b This is a line with slope of -20 DB/decade with a value of 0 DB at omega = 1 rad/sec. Third, there is a pole on the real axis at -B given by the term p15_7c The magnitude of this at low frequencies is:

     p15_8    (omegagtb)
 
and at high frequencies the magnitude is:
p15_9 (omegagtb)



DETERMINATION OF THE PARAMETERS FOR THE CART AND MOTOR


The first is a horizontal line at 0 DB. The second is a line of slope -20 DB/decade which crosses 0 DB at omega = B. We thus expect the Bode plot we obtain experimentally to be the sum of the three individual responses. We expect it to look something like the following:


p16_1


DETERMINATION OF PARAMETERS FOR CART AND MOTOR: SQUARE WAVE INPUT

Square Wave Input:     For omegaltb

Due to static friction inherant in the cart and motor system it is difficult to obtain an accurate Bode plot when pure sinusoidal inputs are used. When the voltage to the motor is near zero the system tends to stick and not move at all. To overcome this problem we will use a square wave input instead of a sine wave input. We now examine the response we expect to get using a square wave input.
A square wave of frequency omega can be expressed as the sum of sine waves as follows:

           p17_1

Since our transfer function for the cart and motor is linear, the output due to a square wave is the sum of the outputs due to each of the sinusoudal terms which make up the square wave.

From our earlier analysis we know the magnitude of the output oscillations due to an input of asinwt is:

           p17_2

For  omegaltb this becomes

           p17_3

For a square wave input expressed as the sum of a series of sinusoidal inputs the magnitude of the output becomes:

  p17_4a p17_4c
 
p17_5c
 
Note:  p17_5d  (Proof)
 
p17_6c

We thus expect the output due to a square wave input to be about 1.233 times the output due to a sine wave input when ( omegaltb)

Expressed in decibels this becomes:

      DB   =   p17_7c
 
  =   p17_8c
 
  =   p17_9c


The output due to a square wave input is about 1.82 DB greater than it would be for a sine wave input when ( omegaltb)


Solution to series 1 + 1/(3^2) + 1/(5^2) + 1/(7^2) + 1/(9^2) + ...
> solution to this series:
>
> 1 + 1/(3^2) + 1/(5^2) + 1/(7^2) + 1/(9^2) + ...
See: http://www.efunda.com/math/seriesofconst/IntRecipSeries.cfm


Or:
Let A = 1 + 1/(3^2) + 1/(5^2) + 1/(7^2) + 1/(9^2) + ...

Using the Riemann's Zeta Function and its 'alternating' form
 Zt(s) = Sum_{n=1}^\infty [ 1/n^s ]
 Za(s) = Sum_{n=1}^\infty [ (-1)^{n-1} / n^s ]
       = [ 1 - 2^{1-s} ] Zt(s)

It is easy to verify that
 A = [ Zt(2) + Za(2) ] / 2
   = (1/8)pi^2

There are around half-dozen proofs that Zeta(2) = pi^2/6.

========================
Or:
Here is a proof using Fourier series:
Take the function fx(x) = -x when -π < x < 0
= 0 when 0 < x < π
Notice f(x) is piecewise continuous with period 2π so this means we can find its corresponding Fourier series.

Solve for the coefficients of the Fourier series, and write out the Fourier series for f(x)

f(x) = π/4 - (2/π)*[cos(x) + (1/9)cos(3x) + (1/25)cos(5x) ...] - [sin(x) - (1/2)sin(2x) + (1/3)sin(3x) ...]

Now just evaluate at x=0

You get

0 = π/4 - 2/π * (1 + 1/9 + 1/25 +...)

so π2/8 = (1 + 1/9 + 1/25 + ...)




Square Wave Input:    For  omegagtb


We now examine what happens when  omegagtb    For  omegagtb our output due to a sine wave input becomes:

          p18_10

For a square wave input expressed as the sum of a series of sinusoidal inputs the output becomes:

      p17_4a       p18_11c
 
p18_12c
 
p18_13c


We thus expect the output due to a square wave input to be about 1.052 times the output due to a sine wave input. Expressed in decibels this becomes:

      DB   =   p18_14c
 
  =   p18_15c
 
  =   p18_16c


The output due to a square wave input is about 0.44 DB greater than it would be for a sine wave input when  omegagtb

We may conclude that the difference in response between a sine wave input and a square wave input is small. Thus, to experimentally determine the unknown parameter B we may use square waves of various frequencies and analyze the results the same as we would for a sine wave input.
BLOCK DIAGRAM


          p19_1


TRANSFER FUNCTIONS


      p11d4       (The Inverted Pendulum)
 
p3_4 (The Cart and Motor System)





We will later use the variables xsubp(t),  xdot1p(t),  theta(t),  thetad1(t) as the state space variables for our state space representation of the system. We can show the Laplace transform of these variables on the block diagram as follows:

          p19_2




Derivation of State Space Representation


Equation of motion for the cart

      p11c6       (Our characteristic equation of motion for the cart and motor system derived earlier.)
p20a2 (Solve for xdot2p)



Equation of motion for the Inverted Pendulum

      p2_9       (Our characteristic equation of motion for the inverted pendulum derived earlier. Now solve for thetad2 )
p20b2
p20b3
p20b4


State Space Representations

      p20c1      
p20a2 (Equation of motion for the cart derived above)
p20c3
p20b4 (Equation of motion for the inverted pendulum derived above)

This can be represented in matrix form as:

p20c5


This has the general form of:

          p20d1



Calculation of Transfer Functions from State Space Equation


For the moment we assume all four of the state space variables are individually observable. Thus:


      p21a1


      This has the form   p21a2   where C = I, the identity matrix


The transfer function matrix is given by:

F(s) = Cp21siaB



Calculate p21sia where I is the identity matrix, and A was the 4x4 matrix derived two pages ago under Derivation of State Space Representation.




p21b1



Calculate F(s) = Cp21siaB where C = I, the identity matrix, and B is the vector derived three pages ago at the bottom under Derivation of State Space Representation.




           F(s) = Cp21siaB


p22_1



p22_2



Recall from earlier the transfer function equations:


p11d4       Transfer function for the cart and motor system.
 
p3_4 Transfer function for the inverted pendulum.


And recall from the Laplace transform that p22_7 = s laplace [f(t)]

Thus p22_4 and p22_5


Our equation F(s) is now:


p22_6
Discretization Of Continuous Time State Space Equations



where

          p23a3

          p23a4


Calculation of G


          p23a3

          p23b2



p23c1




Calculation of inverse Laplace Transforms


We now have the arduous task of calculating the inverse Laplace transform for the terms in the 4x4 matrix on the previous page.

Going down the rows we have

Row 1, Column 1

           pm11_1



Row 1, Column 2


           pm12_2




Row 2, Column 2


           pm22



Calculation of inverse Laplace Transform for row 3 column 2.



         pm32_1


when asub1 = 0 and asub2 = kpb

Substituting these values of asub1 and asub2 into equations N1, N2, N3 on the previous page gives

           pm32_2a

           pm32_2b

           pm32_2c


Substituting these values of nsub1 , nsub2 , and nsub3 into equation L4 on the previous page gives the Laplace transform of row 3 column 2 as

pm32

which reduces to

pm32r



Calculation of inverse Laplace Transform for row 4 column 2.



         pm42_1


when asub1 = 0 and asub2 = ksubpb

Substituting these values of asub1 and asub2 into equations N1, N2, N3 derived two pages ago gives

           pm42_2a

           pm42_2b

           pm42_2c


Substituting these values of nsub1 , nsub2 , and nsub3 into equation L4 derived two pages ago gives the Laplace transform of row 4 column 2 as

pm42

which reduces to

pm42r


Calculation of inverse Laplace Transform for (row,column) (3,3) (3,4)(4,3) (4,4)



A general form of these equations is

           pm33g_4

(We will save nsub1 for later use.)

Extract nsub2 from top and factor the bottom

           pm33g_5

Add zomegan - zomegan to top

           pm33g_6

Split top

           pm33g_7

Multiply second term by pm33g_8j

           pm33g_8

Taking the inverse Laplace transform of the above equation gives

           pm33g_9

EQUATION L5

           pm33g_10


Calculation of inverse Laplace Transform for row 3 column 3.



         pm33_1


when nsub2 = 1 and nsub3 = 2zomegan

Substituting these values of nsub2 and nsub3 into equation L5 derived on the previous page gives the Laplace transform of row 3 column 3 as

           pm33

which slightly reduces to

           pm33s


Calculation of inverse Laplace Transform for row 3 column 4.



         pm34_1


when nsub2 = 0 and nsub3 = 1

Substituting these values of nsub2 and nsub3 into equation L5 gives the Laplace transform of row 3 column 4 as

           pm34

or

           pm34r


Calculation of inverse Laplace Transform for row 4 column 3.



         pm43_1


when nsub2 = 0 and nsub3 = omegansq

Substituting these values of nsub2 and nsub3 into equation L5 gives the Laplace transform of row 4 column 3 as

           pm43

or

           pm43r


Calculation of inverse Laplace Transform for row 4 column 4.



         pm44_1


when nsub2 = 1 and nsub3 = 0

Substituting these values of nsub2 and nsub3 into equation L5 gives the Laplace transform of row 4 column 4 as

           pm44


Calculation of G


           aeq

          p23a3

          p23b2



           draft22

where
(3,2) = pm32r

(4,2) = pm42r

(3,3) = pm33rs

(3,4) = pm34s

(4,3) = pm43s

(4,4) = pm33rs



Calculation of G assuming zeta = 0


If our damping factor for the pendulum zeta is estimated to be small enough that we can set it to zero then our equations simplify to

           draft23

where
(3,2) = pmz32

(4,2) = pmz42r



Calculation of H when zeta, = 0


          p30_1a                       p30_1b



p30_2

rearanging gives

p30_3





p31_1


p31_2

p31_3

p31_4


Discussion of Observable Properties


The two state veriables X and theta are directly observable by measuring the voltage from the corresponding potentiometers and multiplying by appropriate scaling constants. Thus:

      p32a_1a       ;       vsubx = Voltage corresponding to X position       p32a_1c
      p32a_2a       ;       vsubx = Voltage corresponding to X position       p32a_2c


For the moment we will assume the other variables xdot and thetad1 are also directly observable. Thus our output vector Y consisting of observable states of the system becomes:

          p32a_3

which has the form

          p32a_4


Discussion of Feedback


We now decide to use a linear combination of the observable states contained in vector y to determine the input controlling voltage.

           p32b_1

           vsubi = KY
           p32b_3

Our original state space equation was

           x(k+1) = Gx(k) + Hv(k)

           y = Cx(k)

Now observe that our input v(k) is actually vsubi , the input controlling voltage to the servo motor. We also have decided to set vsubi = KY, and we know y = Cx(k) Thus we now have:

           x(k+1) = Gx(k) + HkCx(k)

           x(k+1) = (G + HkC)x(k)

p33_5

Note that in this case HK is the outer product of two vectors producing a square matrix.

If we assume an initial x vector x(0), we then have

p33_6

As k increases, the components of vector x(k) will tend toward zero only if the eigenvalues of the matrix (G + HKC) all have magnitude less than one. It is now our job to chose the four values, p33_7 such that the matrix (G + HKC) is stable.



Procedure to obtain ksubx


Measure voltage vsubl from position potentiometer. Move cart one meter to right and measure voltage vsubr from position potentiometer.

           p34a_1



Procedure to obtain ktheta


Lay pendulum at -90 degrees horizontal with table, pendulum pointing left. Measure vsubl from pendulum angle potentiometer. Turn pendulum 180 degrees ( pi radians) clockwise (positive sense) so pendulum is horizontal with table, pointing right. Measure vsubr from pendulum angle potentiometer.

           p34b_1



BLOCK DIAGRAM OF SYSTEM WITH DISCRETE-TIME FEEDBACK


P35_1



Motor Dead Range


Due to static friction the servo motor has an input voltage threshold below which the motor will not move. The friction is primarily due to the brushes within the servo motor. A graph of the input voltage versus cart velocity might look something like the following:



p36_1

To compensate for this non-linearity we can add to our calculated control voltage the corresponding threshold voltage. If the calculated control voltage is positive we add the threshold voltage to the control voltage to compensate for the deadrange, and if the calculated voltage is negative we add a negative threshold voltage to the control voltage to compensate for the deadrange. A graph of the transfer function would look something like:



p36_2

Note that this function is the inverse of the previous function. It thus cancels out the non-linearity of the motor-cart system due to static friction making the system appear more linear.





BLOCK DIAGRAM OF SYSTEM WITH DISCRETE-TIME FEEDBACK AND MOTOR DEADSPACE CANCELATION



P37_1


PENDULUM SETUP



p38



Parameters for the Cart and Motor


The response of the cart and motor system to square wave inputs of frequencies 0.1 Hz to 15 Hz was recorded on the chart recorder and plotted on a Bode plot. The primary breakpoint between the line of slope -20 Db/decade and the line of -40 Db/decade occurred at 1 Hz (or omega = 6.28 rad.) A secondary breakpoint of unknown origin occurred at 7 Hz.



The cart gain constant ksubp was obtained by analyzing the magnitude of the cart oscillations for low frequencies. Using the equation



           p40_1

the value of ksubc was determined by assuming the cart acceleration xdot2p was zero for low frequencies and then calculating the average velocity xdot1p from the chart recorder.



COMPLETE LIST OF PARAMETERS


The following parameters were determined by experiment:

omegad = 5.41 rad/sec
zeta = 0
omegan = 5.41 rad/sec
ksubp = 3 1/meter
beta = 6.28 rad/sec
ksubc = 1.5 msqvolt
ksubx = 6.75 volts/meter
ktheta = 3.34 volts/meter
T = 0.06 1/sec



CONTROL OF THE PENDULUM



Determination of Feedback Control Parameters

Two computer programs were written to aid in the selection of feedback control values. The first program, ACKERMAN, calculates the control vector K given the desired pole positions using Ackerman's formula. The second program, PENDULUM, does the reverse; it calculates the pole positions given a control vector K.

Control of the Pendulum


Various values of the control vector K were obtained using program ACKERMAN. The final control vector used was X = 3.5, xdot = 7.0, theta = 25.0, thetad1 = 3.0 Using these values the poles were calculated to be at 0.752 plusormi 0.382j, 0.929, 0.982

Using these values the pendulum remained upright with the cart slowly oscillating back and forth covering a distance of about 25 cm. The cart was run for a period of 2 minutes and showed no signs of instability.


Analysis of the Results

In a final effort to justify the endless hours spent in deriving a theoretical analysis of the pendulum system, it is gratifying to be able to sit back and say IT WORKS!

(If you've been paying attention you may have noticed that two of the poles were at 0.929 and 0.982, which is pretty close to being unstable. The reason is I initially made a mistake calculating the control values because I used the full length of the rod for l when l was originally defined as 1/2 the length of the rod. Recalculating using the correct value for l gave the true poles.)



Links For Further Research

There's been further research since this paper was originally written in 1987. My thanks to Reid Sweatman for providing the following comments and links:



Curious side note: this problem is susceptible to fairly easy solution by neural net, and as such, is an example problem in several texts, sometimes slightly modified to balance the stick upright, rather than oscillate it. Of course, this removes a number of exponential terms from the differential solution, but, strangely enough, that doesn't make much difference to the neural net solution, beyond requiring extra feed-forward layers.

I believe a genetic algorithm might also be capable of handling the problem, so long as it was updated often enough to keep the stick reasonably within the solution space. However, I've never seen that done, or tried it myself. I do neural nets more frequently than genetic algorithms, but have a friend who has done several exceedingly complex game AIs (essentially very complex routing schedulers) who maintains it's possible. I suppose of equal interest is the frequency with which AIs are simulated by physical systems, like constrained springs.

I've got two favorite books on neural nets:

For a web site, you can't do much better for a launch point than Steve Woodcock’s site,. Lots of good stuff there. Lots of really funny things (check the "You know your game is in trouble when..." page), and links to just about everything in AI. Although the site is dedicated to AI in computer games, Steve works for Raytheon doing, last I talked to him, missle guidance systems. Before his brief sojourn into gaming, he wrote AI code for Star Wars (not the movie one). I think he still co-chairs the AI roundtable discussions at the annual Computer Game Developers' conference in San Jose, as well as a more academic one just up the road about two weeks later. Those roundtables are the thing I always like best about the GDC, although I haven't been since 2000. Hey, where else can you get a few hours of what's essentially free consulting with guys like Steve and Neil Kirby from Lucent on your game AI? Hashed out a messy flocking/steering infantry AI with Steve, then the same for training a neural net flight simulation pilot later that night with Neil. Keeps the old brain cells in motion.