A Flaw in Lagrangian Mechanics for a rotary inverted pendulum system? (Simulation vs. Measured)

[Remusco]

TL;DR Summary

I built a rotary inverted pendulum system, and am getting different natural frequencies in the simulation vs. in reality.

I have built this inverted pendulum system (powered by a stepper motor).

I define theta to be the pendulum angle 0 degrees in the vertical position, phi is the motor arm angle.
When I let the physical pendulum go at theta = 90 degrees, with the stepper motor turned off, I get:

Now if I also plot the simulated result on top of this (with equations derived from Lagrangian Mechanics, I get:

However, the frequencies do line up when I use basic Newtonian mechanics. Why is there a difference?

This is the equation I derived using Lagrangian mechanics:

To find b_theta I simply iterated until both plots aligned.

As another weird thing. If I just use I_p in the denominator, the pendulum oscillates correctly (with the correct frequency). However, then the pendulum performs horribly when any control torque is applied. So is there a middle ground here that is correct? Do I use I_p or 4*I_p + m_p*L_p**2???? Either way does not make sense. My derivation is in the Colab link below.

https://colab.research.google.com/drive/1_uAv3Jj6wQPyN1OhDRKam-ZLRShalRYy?usp=sharing

 

[PeroK]

Lagrangian mechanics is equivalent to Newtonian mechanics. If you get different answers, you've made an error somewhere.

 

[renormalize]

TL;DR Summary: I built a rotary inverted pendulum system, and am getting different natural frequencies in the simulation vs. in reality.

Your rotary inverted pendulum is also known as a Furuta pendulum. An analysis of this system, including the inertia tensors, the Lagrangian, the non-linear equations of motion and numerical modeling results, is provided in this paper: On the Dynamics of the Furuta Pendulum

Try comparing your results to those in the paper.

 

[Remusco]

Your rotary inverted pendulum is also known as a Furuta pendulum. An analysis of this system, including the inertia tensors, the Lagrangian, the non-linear equations of motion and numerical modeling results, is provided in this paper: On the Dynamics of the Furuta Pendulum
View attachment 359731
Try comparing your results to those in the paper.

Thanks for sending that to me. Great resource. However, I used the exact results in the paper and this is what i get. Still the same thing is happening:

import sympy as sp

# === Step 1: Define all symbolic variables ===
theta1, theta2 = sp.symbols('theta1 theta2')
theta1_dot, theta2_dot = sp.symbols('theta1_dot theta2_dot')
theta1_ddot, theta2_ddot = sp.symbols('theta1_ddot theta2_ddot')
tau1, tau2 = sp.symbols('tau1 tau2')
J0_hat, J2_hat = sp.symbols('J0_hat J2_hat')
J0, J1, J2 = sp.symbols('J0, J1, J2')
g = sp.symbols('g')

# Constants
m1, m2, L1, L2 = sp.symbols('m1 m2 L1 L2')
b1, b2 = sp.symbols('b1 b2')

# === Step 2: Numerator Expression ===
T = sp.Matrix([theta1_dot, theta2_dot, theta1_dot * theta2_dot, theta1_dot**2, theta2_dot**2])

A = sp.Matrix([
    [m2*L1*(0.5*L2)*sp.cos(theta2)*b1],
    [-b2 * (J0_hat + J2_hat * sp.sin(theta2)**2)],
    [m2*L1*(L2*0.5)*J2_hat * sp.cos(theta2) * sp.sin(2*theta2)],
    [-(1/2)*sp.sin(2*theta2)*(J0_hat*J2_hat + (J2_hat**2) * sp.sin(theta2)**2)],
    [-(1/2)*m2**2 * L1**2 * (L2/2)**2 * sp.sin(2*theta2)]
])

B = sp.Matrix([
    [-m2 * L1 * (L2*0.5) * sp.cos(theta2)],
    [J0_hat + J2_hat * sp.sin(theta2)**2],
    [-m2 * (L2*0.5) * sp.sin(theta2) * (J0_hat + J2_hat * sp.sin(theta2)**2)]
])

U = sp.Matrix([tau1, tau2, g])

numerator = (A.T * T)[0] + (B.T * U)[0]

# === Step 3: Denominator Expression ===
denominator = (J0_hat*J2_hat + J2_hat**2 * sp.sin(theta2)**2) - m2**2 * L1**2 * (L2/2)**2 * sp.cos(theta2)**2

# === Step 4: Final expression ===
theta2_ddot_expr = numerator / denominator

theta2_ddot_expr = theta2_ddot_expr.subs(J0_hat, J1 + m1*(L1/2)**2 + m2*L1**2)
theta2_ddot_expr = theta2_ddot_expr.subs(J2_hat, J2 + m2*(L2/2)**2)


# === Step 5: Lambdify ===
vars = (theta1, theta2, theta1_dot, theta2_dot, tau1, tau2, g, m1,
        m2, L1, L2, b1, b2, J1, J2)

thetaddot_func = sp.lambdify(vars, theta2_ddot_expr, 'numpy')

theta2_ddot_expr

When I use a Newtonian model I get:

 

[renormalize]

Thanks for sending that to me. Great resource. However, I used the exact results in the paper and this is what i get. Still the same thing is happening:

Can you explain in more detail how you determined that "the same thing is happening"? I ask because in your link you state:
"The rotary inverted pendulum system is a 1 DOF system with the motor actuator. It will have one Lagrange Equation."
and indeed you find one Lagrange equation-of-motion for your angle $\theta$ and your other angle $\phi$ apparently plays no role in your dynamics.
In contrast, the paper I referenced clearly describes a 2 DOF system and arrives at two coupled Lagrange equations-of-motion for the coordinates $\theta_1,\theta_2$ (your $\phi,\theta\,$):

Why is there no dynamics for $\phi$ in your approach and how are you comparing the results?

 

[Remusco]

Can you explain in more detail how you determined that "the same thing is happening"? I ask because in your link you state:
"The rotary inverted pendulum system is a 1 DOF system with the motor actuator. It will have one Lagrange Equation."
and indeed you find one Lagrange equation-of-motion for your angle $\theta$ and your other angle $\phi$ apparently plays no role in your dynamics.
In contrast, the paper I referenced clearly describes a 2 DOF system and arrives at two coupled Lagrange equations-of-motion for the coordinates $\theta_1,\theta_2$ (your $\phi,\theta\,$):
View attachment 359750
Why is there no dynamics for $\phi$ in your approach and how are you comparing the results?

I'm directly controlling phi with my stepper motor. The position, velocity, acceleration. The dynamics of phi are whatever I want them to be, and I'm assuming that the stepper is rigid and does not slow down under changes in torque, which is a valid assumption. This reduces the DOF by one; therefore I can just sub in the details for phi. When I put phi_dot =0, and let the pendulum oscillate, I get a slower frequency than what should be in reality.

 

[renormalize]

This reduces the DOF by one; therefore I can just sub in the details for phi. When I put phi_dot =0, and let the pendulum oscillate, I get a slower frequency than what should be in reality.

Thanks for the explanation. Can you post in LaTeX the Newton and Lagrange equations you get for $\theta$ to help us understand why they predict different frequencies?

 

[Remusco]

Thanks for the explanation. Can you post in LaTeX the Newton and Lagrange equations you get for $\theta$ to help us understand why they predict different frequencies?

Here are the Newton's equations:
$$ -T \sin{\left(\theta \right)} - b \theta_{dot} \cos{\left(\theta \right)} - m_{p} x_{ddot} - \frac{1}{2} L m_{p} \cos{\left(\theta \right)} \theta_{ddot} + \frac{1}{2} L m_{p} \theta_{dot}^{2} \sin{\left(\theta \right)} = 0$$

$$-T \cos{\left(\theta \right)} - g m_{p} + b \theta_{dot} \sin{\left(\theta \right)} + \frac{1}{2} L m_{p} \theta_{ddot} \sin{\left(\theta \right)} + \frac{1}{2} L m_{p} \theta_{dot}^{2} \cos{\left(\theta \right)} = 0$$
$$ x_{ddot} = L_a *\phi_{ddot}$$


Note that these can be reduced to one equation if we take phi to be driven by an actuator.
I ended up using Newton equations for my simulation since I couldn't get the Lagrange equation to work.

Here is the Lagrange equation:

$$\ddot{\theta} = \frac{0.5 \left(-4.0 L_a L_p m_p \dot{\phi} \cos(\theta) + L_p^2 m_p \dot{\phi}^2 \sin(2\theta) + 4.0 L_p g m_p \sin(\theta) - 8.0 b_\theta \dot{\theta} \right)}{4.0 I_p + L_p^2 m_p}$$

Below is my result with the Newton derivation. Although the inertia seems underestimated, I'm happy with the result overall.

 

[renormalize]

Here are the Newton's equations:$$ -T \sin{\left(\theta \right)} - b \theta_{dot} \cos{\left(\theta \right)} - m_{p} x_{ddot} - \frac{1}{2} L m_{p} \cos{\left(\theta \right)} \theta_{ddot} + \frac{1}{2} L m_{p} \theta_{dot}^{2} \sin{\left(\theta \right)} = 0$$$$-T \cos{\left(\theta \right)} - g m_{p} + b \theta_{dot} \sin{\left(\theta \right)} + \frac{1}{2} L m_{p} \theta_{ddot} \sin{\left(\theta \right)} + \frac{1}{2} L m_{p} \theta_{dot}^{2} \cos{\left(\theta \right)} = 0$$

I've rewritten these in the more conventional form:$$-T\sin\theta-b\,\dot{\theta}\cos\theta-m_{p}\ddot{x}-\frac{1}{2}Lm_{p}\cos\theta\,\ddot{\theta}+\frac{1}{2}Lm_{p}\dot{\theta}^{2}\sin\theta=0$$$$-T\cos\theta-gm_{p}+b\,\dot{\theta}\sin\theta+\frac{1}{2}Lm_{p}\ddot{\theta}\sin\theta+\frac{1}{2}Lm_{p}\dot{\theta}^{2}\cos\theta=0$$Can you verify that I've gotten your equations correct? And what is the quantity $\ddot{x}$ that appears in the first one?

Note that these can be reduced to one equation if we take phi to be driven by an actuator.

I don't see $\phi,\dot{\phi},\ddot{\phi}$ anywhere in your equations. Or are they somehow in $T\,$? (And is that a torque?)

Here is the Lagrange equation:$$\ddot{\theta} = \frac{0.5 \left(-4.0 L_a L_p m_p \dot{\phi} \cos(\theta) + L_p^2 m_p \dot{\phi}^2 \sin(2\theta) + 4.0 L_p g m_p \sin(\theta) - 8.0 b_\theta \dot{\theta} \right)}{4.0 I_p + L_p^2 m_p}$$

Can you confirm that the decimal numbers appearing here represent the exact values $1/2,4,8\,$?

Below is my result with the Newton derivation. Although the inertia seems underestimated, I'm happy with the result overall.
View attachment 359755

Your "Simulation" and "Physical Model" labels confuse me. Presumably one is the solution to your Newtonian equation(s). Which one is it? And how exactly are you generating the other curve that you compare against? Frankly, to me neither of the plots exhibit particularly good agreement.

 

[Remusco]

I've rewritten these in the more conventional form:$$-T\sin\theta-b\,\dot{\theta}\cos\theta-m_{p}\ddot{x}-\frac{1}{2}Lm_{p}\cos\theta\,\ddot{\theta}+\frac{1}{2}Lm_{p}\dot{\theta}^{2}\sin\theta=0$$$$-T\cos\theta-gm_{p}+b\,\dot{\theta}\sin\theta+\frac{1}{2}Lm_{p}\ddot{\theta}\sin\theta+\frac{1}{2}Lm_{p}\dot{\theta}^{2}\cos\theta=0$$Can you verify that I've gotten your equations correct? And what is the quantity $\ddot{x}$ that appears in the first one?

I don't see $\phi,\dot{\phi},\ddot{\phi}$ anywhere in your equations. Or are they somehow in $T\,$? (And is that a torque?)

Can you confirm that the decimal numbers appearing here represent the exact values $1/2,4,8\,$?

Your "Simulation" and "Physical Model" labels confuse me. Presumably one is the solution to your Newtonian equation(s). Which one is it? And how exactly are you generating the other curve that you compare against? Frankly, to me neither of the plots exhibit particularly good agreement.

The physical model result is from my physical rotary inverted pendulum. It is the data I collected with my Arduino. The simulated as shown in the plot is the result from the Newtonian equations. Here I now solve the Newtonian equations for theta_ddot:

 

[renormalize]

The simulated as shown in the plot is the result from the Newtonian equations. Here I now solve the Newtonian equations for theta_ddot:
View attachment 359781

I'm forced again to rewrite this in LaTeX:$$\ddot{\theta}=\frac{2\left(-b\,\dot{\theta}+gm_{p}\sin\theta-L_{a}\ddot{\phi}\cos\theta\right)}{L_{p}m_{p}}$$I fail to understand how this equation is obtained from your 2 Newton's eqs. that I rewrote in post #9, especially since $T$ and $\ddot{x}$ in those are nowhere defined and $\phi$ is evidently absent entirely. I'm sorry, your presentation of the problem is simply too scattershot and incoherent for me to follow, so I'm dropping out of this conversation.
Good luck.

 

[Remusco]

I'm forced again to rewrite this in LaTeX:$$\ddot{\theta}=\frac{2\left(-b\,\dot{\theta}+gm_{p}\sin\theta-L_{a}\ddot{\phi}\cos\theta\right)}{L_{p}m_{p}}$$I fail to understand how this equation is obtained from your 2 Newton's eqs. that I rewrote in post #9, especially since $T$ and $\ddot{x}$ in those are nowhere defined and $\phi$ is evidently absent entirely. I'm sorry, your presentation of the problem is simply too scattershot and incoherent for me to follow, so I'm dropping out of this conversation.
Good luck.

All I'm trying to do is replicate this YouTube tutorial, but for a uniform pendulum rod, not a point mass.
https://www.google.com/search?q=inv...ate=ive&vld=cid:93075d53,vid:5qJY-ZaKSic,st:0