How to develop a physically correct simulation of an inverted pendulum

Start with a picture of the system we want to model.

Finding the Equations of Motion

We want to develop a set of equations that describe the dynamics of this pendulum. These equations are known as the equations of motion. To get these equations I prefer to use the Lagrangian[1] method because it is scalable to much more complex simulations.

The Lagrangian is L = T- U where T is the kinetic energy of the system and U is the potential energy.

T is the sum of kinetic energies of m₁ and m₂

T = 1/2 m₁ v₁² + 1/2 m₂ v₂²

where v₁ and v₂ are the velocities of m₁ and m₂ respectively.

U = -m₁ g l cos(Φ)

Giving us the Lagrangian

L = 1/2 m₁ v₁² + 1/2 m₂ v₂² + m₁ g l cos(Φ)

We have to get L in terms of l, x and Φ instead of v₁ and v₂. I will use Mathematica to work through the equations, as it is a very nice tool for such work. The details are here, and give the resulting Lagrangian as

L = l cos(Φ(t)) ( g m₁ + m₂ x'(t) Φ'(t) ) + 1/2 (m₁x'(t)² + m₂ ( x'(t)² + l² Φ'(t)² ) )

The Lagrangian method is nice because given this equation we can now find the equations of motion by solving the following equation.[1],[2]

Where q_k are the generalized coordinates of the system and Q_k the generalized force, in our case x, Φ, F_x, and F_Φ. From this we find two coupled differential equations of motion.

x''(t) = ( F_x(t) + l m₂ ( sin(Φ(t)) Φ'(t)² - cos(Φ(t)) Φ''(t) )) / ( m₁ + m₂ )

Φ''(t) = ( F_ø(t) - l ( g m₁ sin(Φ(t)) + m₂ cos(Φ(t)) x''(t) ) ) / l² m₂

Integrating the Equations of Motion

Now that we have equations describing the dynamic motion of our pendulum we need to solve them. As our equations are differential equations, solving involves integration of the equations. In general anything but the simplest differential equations can not be integrated in closed form. That is to say there is no simple solution involving known functions like sine or cosine. Since our intention is to come up with a computer simulation we will use a numerical method of integration. Runge-Kutta[3] is the numerical method that I like to use. Runge-Kutta works with first order differential equations, and we will use the 4th order Runge-Kutta which is computed as follows.

Given dy/dx = f(x,y), we want to find y(x).

Suppose y = y0 at x = x0, we choose a small interval h and calculate

k1 = h f(x0, y0)
k2 = h f(x0 + h/2, y0 + k1/2)
k3 = h f(x0 + h/2, y0 + k2/2)
k4 = h f(x0 + h, y0 + k3)

Then at x = x0 + h , y = y0 + 1/6(k1 + 2 k2 + 2 k3 + k4).

Since we have second order equations, we first need to put them in first order form.
Noting that x''(t) = v'(t) and Φ''(t) = ω'(t) we substitute and solve our equations of motion for v'(t), x'(t), Φ'(t) and ω'(t) giving us four single order equations.

x'(t) = v(t)

Φ'(t) = ω(t)

v'(t) = (2(l F_x(t)-cos(Φ(t))F_Φ(t)+l sin(Φ(t))(gml cos(Φ(t)) +lm₂ω(t)²))) /
(l (2m₁+m₂-m₂ cos(2Φ(t))) )

Now that we have our equations in the right form we can write the software.

Model View Controller

We will write a Java applet using the model view controller pattern, or MVC[4]. The MVC breaks up the problem into the three named classes.

Model

The model is responsible for simulating the behavior of the system. This is where we use the derived equations of motion and where our Runge-Kutta method is implemented. The model needs to keep track of the state of the system, to iterate the system through small time steps, and allow us to get the results.

We will start with a Java class IPModel.

• reset() - puts the model in an initial known state.
• toString() - return a String representation of the model state.
• update(double deltaTime) - propagates the model through deltaTime using the Runge-Kutta method.
• getX() - return current value of x
• getPhi() - return current value of Φ
• getL() - return current value of l

A more complete implementation would provide get() and set() methods for all the state variables but we do not need them for this demo and they are left out for clarity.

You can compile and run IPModel.java by itself. The main() method will run the model through 10 seconds of time and output the model state to the console.

View

The view needs to display the block, the pendulum mass, and the pendulum rod. We will subclass javax.swing.JPanel as our IPView class.

• setModel(IPModel) - associate the view with a specific model instance.
• modelStateChanged() - tell the view to update itself to reflect a new model state.

The protected method paintComponent() handles the actual drawing. This method is not called directly, but indirectly via the this.repaint() in modelStateChanged().

A note about the Netbeans IDE (Integrated Development Environment).
I use the Netbeans IDE for Java development which includes a GUI (Graphical User Interface) editor for developing custom views. Looking at the code for IPView and IPApplet you may notice comments warning about modification, specifically to initComponents(). This code is generated by the GUI editor and not hand written. For our purposes it can be ignored, but it is worth looking at if you want to fully understand the construction of the Applet.

Controller

Our controller is fairly simple. At a regular interval of 1/10 second the controller updates the model and the view. Therefore we do not have a separate controller class but integrate the functionality of the controller into the IPApplet class. The method init() creates an instance of IPModel and associates it with our IPView and finally starts a timer which calls actionPerformed() periodically. In turn actionPerfomed() updates the model and notifies the view that the model state has changed.

Here is the result.

IPApplet

Where to go from here

This demo is complete in the sense that we have an accurate representation of the physics of the system, but it is limited in that the parameters, m₁, m₂ , and l are fixed. Separating the controller from IPApplet and extending it to manipulate these model parameters would be interesting. IPView should be modified in this case to better represent these values. In this case consideration needs to be taken so that the length l is correctly proportional to the distance x moved by m₁.

Furthermore, the model allows for F_x, and F_Φ to be non zero. One could implement a controller that provides feedback to F_x such that the pendulum is balanced in the vertical position against random fluctuations applied to F_Φ.

Comments

I am interested in any difficulties encountered in using this information.
Please contact me with any comments, bugs, or suggestions.