Error in total energy of the cube experiments data/tosses_processed

[sunilgora]

Hi, I am Sunil from IIT Kanpur, and I am working on contact dynamics to simulate cube toss experiments. The total energy, i.e., the sum of kinetic and potential energy of the cube toss data stored in data/tosses_processed after impact, comes out to be greater than the initial total energy of the cube for some of the experiment data e.g., 0.pt, 7.pt, 12.pt etc. as shown below.

Please check the Python code below to reproduce these plots:

import torch,os
import numpy as np
import matplotlib.pyplot as plt

dt = 1/148 #148 Hz

#Cube parameters (Acosta et al 2022)
m = 0.37 # mass of the cube in kg
I = 0.0081 #mass moment of inertia
width = 0.1 #10x10x10cm cube
mu=0.18 #friction
e=0.125 #restitution
blockWidthScale =0.0524 # scale factor used in dataset

#Get current directory
base_dir = os.path.dirname(os.path.abspath(__file__))

#Load trajectory data .pt from GitHub/contact-nets/data/tosses_processed/file.pt
data_dir=os.path.join(base_dir, "GitHub/contact-nets/data/tosses_processed")

for i in range(0, 20): # files are named 0.pt, 1.pt, ..., 570.pt
    data_path = os.path.join(data_dir, f"{i}.pt")

    #load .pt file data using torch
    data = torch.load(data_path)
    print(f"Loaded data with shape {data.shape}")

    #Remove batch dimension if present
    data = data.squeeze(0)  # Now shape [111, 19]

    #Assume features: 0=x, 1=y, 2=z, 3=q0, 4=q1, 5=q2, 6=q3, 7=dx, 8=dy, 9=dz, 10=wx, 11=wy, 12=wz, 13-18=zeros
    #Time is not stored,
    #frequency is 148 Hz, so time step is 1/148 seconds
    t = dt*np.arange(data.shape[0])

    #plot total energy
    KE= 0.5 * m * blockWidthScale**2*(data[:, 7]**2 + data[:, 8]**2 + data[:, 9]**2) + 0.5 * I * (data[:, 10]**2 + data[:, 11]**2 + data[:, 12]**2)
    PE= m * 9.81 * blockWidthScale*data[:, 2]  # m*g*h
    TE= KE + PE
    plt.figure()
    plt.plot(t, TE, marker='o', label='KE+PE')
    plt.xlabel('Time (s)')
    plt.ylabel('Energy (J)')
    #plt.title('Total Energy during Toss')
    plt.legend()
    plt.grid()

plt.show()

The total energy after impact can not be greater than the initial energy of the cube. Could you check and let me know if there are any issues with the above Python code or the experimental data?

[mposa]

This is a really nice experiment. My guess is that the issue is with the moment of inertia that you used. Moments of inertia are very hard to identify, and these trajectories are so short that the value identified in the Acosta paper is likely not physically accurate (even if it produces accurate trajectories in simulation).

From a back of the envelope calculation, I suspect that the value you used for I is too large by a factor of 2x-5x. Try tweaking it, and see how that effects your energy plot.

[sunilgora]

Yes indeed, the mass moment of inertia is too large, and dividing it by $5$ removes this error in total energy. In Acosta's paper, it is mentioned that 'The cube was weighed to determine the mass, with the inertia determined from the measured mass and geometry', but the given mass moment of inertia value of $0.0081 kg\text{-}m^2$ is much larger than the calculated value of $0.00061 kg\text{-}m^2$ ($I=\frac{1}{6}m \times \text{l}^2$) for $m=0.37 kg$ and $l=0.1 m$. Also, if $I=0.0081 kg\text{-}m^2$ is used directly in MuJoCo/Drake/Bullet simulator for parameter identification, the simulator will not generate the trajectories that violate energy conservation. Wouldn't it give the incorrect trajectory errors for the cube simulation?

[mposa]

On second thought, my guess at the inertia value may have been off (and the number Brian used does feel reasonable). The cube is approximately hollow, so the formula you gave isn't quite right. All that said, obviously energy was conserved in the physical experiments.

• I don't see any immediate error in your code snippet, though perhaps there's something off with how the unit conversions are done (though it looks right to me).
• Moment of inertia is only weakly observable in this type of data, and it's possible that a bad moment of inertia will cause a negligible effect on trajectory performance.

[sunilgora]

The mass moment of inertia estimated from the conservation of energy for each trajectory before impact (taking initial $5$ timesteps data for each trajectory) is also around $0.0033 kg\text{-}m^2$, which is closer to the theoretical value of $0.001 kg\text{-}m^2$ for a hollow cube than the value of $0.0081 kg\text{-}m^2$ used in the paper.
Please check this Python code:

import torch,os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

dt = 1/148 #148 Hz

#Cube parameters (Acosta et al 2022)
width = 0.1 #10x10x10cm cube
m = 0.37 # mass of the cube in kg
I_cube = 5/18*m*width**2 #0.0081/5 #1/6*m*width**2 #0.0081 #mass moment of inertia
print("Theoretical mass moment of inertia:", I_cube)
mu=0.18 #friction
e=0.125 #restitution
blockWidthScale =0.0524 # scale factor used in dataset

# Get current directory
base_dir = os.path.dirname(os.path.abspath(__file__))

# Load trajectory data .pt from GitHub/contact-nets/data/tosses_processed/file.pt
data_dir=os.path.join(base_dir, "GitHub/contact-nets/data/tosses_processed")

n_expdata =570 #number of trajectories in the dataset

COMz_data = np.zeros((n_expdata, 5)) #for sys identification, we want to store the height data for each trajectory in a list of arrays
twist_data = np.zeros((n_expdata, 5, 6)) #for sys identification, we want to store the twist (velocity and angular velocity) data for each trajectory in a list
cube_data=[]
for i in range(0, n_expdata): # Assuming files are named 0.pt, 1.pt, ..., 570.pt
    data_path = os.path.join(data_dir, f"{i}.pt")

    #load .pt file in cube_data using torch
    cube_data.append({'exp{i}': torch.load(data_path).squeeze(0)}) #squeeze to remove batch dimension, now shape [111, 19]
    # print(f"Loaded data with shape {data.shape}")

    # Assume features: 0=x, 1=y, 2=z, 3=q0, 4=q1, 5=q2, 6=q3, 7=dx, 8=dy, 9=dz, 10=wx, 11=wy, 12=wz, 13-18=zeros
    # Time is not stored,
    # frequency is 148 Hz, so time step is 1/148 seconds
    t = dt*np.arange(cube_data[i]['exp{i}'].shape[0])

    #Clip data for energy conservation before contact
    data_en0 = cube_data[i]['exp{i}'][:5, :]

    COMz_data[i, :cube_data[i]['exp{i}'].shape[0]] = blockWidthScale*data_en0[:, 2].numpy() #store height data for sys identification
    twist_data[i, :cube_data[i]['exp{i}'].shape[0], 0:3] = blockWidthScale*data_en0[:, 7:10].numpy() #store twist data for sys identification
    twist_data[i, :cube_data[i]['exp{i}'].shape[0], 3:6] = data_en0[:, 10:13].numpy() #store twist data for sys identification

#Find mass moment of inertia from conservation of total energy data

def my_fun(I):
    #Find TE for each traj and minimize the variance of TE across time for each trajectory, averaged across trajectories
    TE_var = []
    for i in range(n_expdata):
        TE_traj=np.zeros(twist_data[i].shape[0])
        for j in range(len(twist_data[i])):
            KE= 0.5 * m * (twist_data[i, j, 0]**2 + twist_data[i, j, 1]**2 + twist_data[i, j, 2]**2) + 0.5 * I * (twist_data[i, j, 3]**2 + twist_data[i, j, 4]**2 + twist_data[i, j, 5]**2)
            PE= m * 9.81 * COMz_data[i, j]  # m*g*h
            TE_traj[j]= KE[0] + PE
        TE_var.append(np.var(TE_traj))
    return np.mean(TE_var)

res = minimize(my_fun, x0=[0.0])
print("Estimated mass moment of inertia:", res.x[0])
I = res.x[0]