Predict trajectory of a bouncing ball

Key Points:

• I have a default ball trajectory generated using some code(provided below). Lets name this trajectory Initial Trajectory.
• Next I have an actual ball whose trajectory I need to estimate.
• Now I need to improve the Initial Trajectory iteratively (to better represent actual ball trajectory) based on the coordinates of the actual ball as it moves through the environment.
• For example lets say after 5 units of time, I have 5 positions of where the actual ball has been through. Now I need to uses these 5 points and the Initial Trajectory to predict the future trajectory of the actual ball. I am expecting the actual trajectory to improve with time as more positions of the actual ball comes in.

Initial Trajectory Visualization:

Initial Trajectory Code:

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81  # gravity (m/s^2)
e = 0.8  # coefficient of restitution
theta = np.radians(45)  # launch angle (radians)
v0 = 10  # initial velocity (m/s)
x0, y0 = 0, 2  # initial position (m)
vx0 = v0 * np.cos(theta)  # initial horizontal velocity (m/s)
vy0 = v0 * np.sin(theta)  # initial vertical velocity (m/s)
t_step = 0.01  # time step for simulation (s)
t_max = 5  # max time for simulation (s)

# Initialize lists to store trajectory points
x_vals = [x0]
y_vals = [y0]

# Simulation loop
t = 0
x, y = x0, y0
vx, vy = vx0, vy0

while t < t_max:
    t += t_step
    x += vx * t_step
    y += vy * t_step - 0.5 * g * t_step ** 2
    vy -= g * t_step

    if y <= 0:  # Ball hits the ground
        y = 0
        vy = -e * vy

    x_vals.append(x)
    y_vals.append(y)

    if len(x_vals) > 1 and x_vals[-1] == x_vals[-2] and y_vals[-1] == y_vals[-2]:
        break  # Stop simulation if ball stops bouncing

# Plotting the trajectory
plt.figure(figsize=(10, 5))
plt.plot(x_vals, y_vals, label="Trajectory of the Ball")
plt.xlabel("Horizontal Distance (m)")
plt.ylabel("Vertical Distance (m)")
plt.title("Trajectory of a Bouncing Ball")
plt.legend()
plt.grid(True)
plt.show()

QUESTION:

• How can I combine a default bouncing ball trajectory and some coordinates of actual ball trajectory to predict the future trajectory of the actual ball.

NOTE:

• I am going with this method of trajectory approximation to better represent real world ball trajectory without having to fine tune complex parameters like energy loss on ball bounce, friction, type of floor, ball elasticity etc.

EDIT:

Actual Ball Trajectory Data:

X: [7.410000e-03 9.591000e-02 2.844100e-01 5.729100e-01 9.614100e-01
 1.449910e+00 2.038410e+00 2.726910e+00 3.373700e+00 4.040770e+00
 4.800040e+00 5.577610e+00 6.355180e+00 7.132750e+00 7.910320e+00
 8.687900e+00 9.465470e+00 1.020976e+01 1.092333e+01 1.163690e+01
 1.235047e+01 1.306404e+01 1.377762e+01 1.449119e+01]

Y: [2.991964 2.903274 2.716584 2.431894 2.049204 1.568514 0.989824 0.313134
 0.311512 0.646741 0.88397  1.0232   1.064429 1.007658 0.852887 0.600116
 0.249345 0.232557 0.516523 0.702488 0.790454 0.78042  0.672385 0.466351]

Graph(Actual Ball Trajectory):

• Note the ball being dropped from a certain height with some horizontal velocity.
• The graph is in 3D but I have simplified the data points to be 2D only (one of the axis is always zero) to keep the problem simpler.

What you are trying to do is known as system identification in the academic literature. Sys Id is a method used to identify the parameters of dynamical systems models from trajectory data.

Once you have a good model fitted to the data, you can then predict the future trajectory of the ball (see also @el_pazzu's answer for using a Kalman filter although this might be tricky with your non-linear system).

Your system is non-linear due to the bouncing behaviour so you will have to set up a non-linear optimization problem and try to solve it for the parameters of your model.

The simplest method of system identification is the single-shooting method which involves simulating an entire trajectory and then comparing it to the outputs from the true system.

This is the prediction error method, which usually involves minimizing the squared-errors between the output trajectory predicted by the model and the data.

Hope that helps. There are many resources online to do this. There are some curated tutorials here:

• https://resourcium.org/topic/system-identification

Also, MATLAB software has a comprehensive set of tools to help you do it.

But I would try solving it using scipy.optimize first. If that doesn't work you may need to consider more advanced non-linear optimization algorithms (see for example, GEKKO, CasADi, Pyomo

Hope that helps. Good luck, it's a very interesting problem...

Have you tried using a Kalman filter:

import numpy as np
import matplotlib.pyplot as plt

g = 9.81
e = 0.8
t_step = 0.01

# Initial actual ball trajectory data
actual_x = np.array([7.410000e-03, 9.591000e-02, 2.844100e-01, 5.729100e-01, 9.614100e-01,
                     1.449910e+00, 2.038410e+00, 2.726910e+00, 3.373700e+00, 4.040770e+00,
                     4.800040e+00, 5.577610e+00, 6.355180e+00, 7.132750e+00, 7.910320e+00,
                     8.687900e+00, 9.465470e+00, 1.020976e+01, 1.092333e+01, 1.163690e+01,
                     1.235047e+01, 1.306404e+01, 1.377762e+01, 1.449119e+01])
actual_y = np.array([2.991964, 2.903274, 2.716584, 2.431894, 2.049204, 1.568514, 0.989824, 0.313134,
                     0.311512, 0.646741, 0.88397, 1.0232, 1.064429, 1.007658, 0.852887, 0.600116,
                     0.249345, 0.232557, 0.516523, 0.702488, 0.790454, 0.78042, 0.672385, 0.466351])

# Initial state (position & velocity)
initial_state = np.array([0, 2, 10 * np.cos(np.radians(45)), 10 * np.sin(np.radians(45))])  # [x, y, vx, vy]

# State transition matrix
dt = t_step
F = np.array([[1, 0, dt, 0],
              [0, 1, 0, dt],
              [0, 0, 1, 0],
              [0, 0, 0, 1]])

# Measurement matrix
H = np.array([[1, 0, 0, 0],
              [0, 1, 0, 0]])

# Process noise covariance
Q = np.eye(4) * 0.001

# Measure noise covariance
R = np.eye(2) * 0.01

# Initial covariance matrix
P = np.eye(4)

def predict(state, P):
    state_pred = F @ state
    P_pred = F @ P @ F.T + Q
    return state_pred, P_pred

def update(state_pred, P_pred, measurement):
    y = measurement - H @ state_pred
    S = H @ P_pred @ H.T + R
    K = P_pred @ H.T @ np.linalg.inv(S)
    state_upd = state_pred + K @ y
    P_upd = (np.eye(4) - K @ H) @ P_pred
    return state_upd, P_upd

# Initialize state & covariance
state = initial_state
trajectory = [state[:2]]
time = 0

# Do Kalman filtering with actual measurements
for i in range(len(actual_x)):
    # Predict
    state_pred, P_pred = predict(state, P)

    # Measure
    measurement = np.array([actual_x[i], actual_y[i]])

    # Update
    state, P = update(state_pred, P_pred, measurement)

    # Save
    trajectory.append(state[:2])

# Convert trajectory to np array for plotting
trajectory = np.array(trajectory)

# Plot trajectory
plt.figure(figsize=(10, 5))
plt.plot(trajectory[:, 0], trajectory[:, 1], label="Predicted trajectory", color='r')
plt.scatter(actual_x, actual_y, label="Actual trajectory", color='b')
plt.xlabel("Horizontal distance")
plt.ylabel("Vertical distance")
plt.title("Trajectory of bouncing ball with Kalman filter")
plt.legend()
plt.grid(True)
plt.show()

This code uses the Kalman filter to refine the trajectory prediction iteratively as more actual data points become available.

I obtained this plot:

You can calculate the discrete second derivative of Y to identify parabolic sections of motion.

from numpy import *
y=array([2.991964, 2.903274, 2.716584, 2.431894, 2.049204, 1.568514, 0.989824, 0.313134, 0.311512, 0.646741, 0.88397, 1.0232, 1.064429, 1.007658, 0.852887, 0.600116, 0.249345, 0.232557, 0.516523, 0.702488, 0.790454, 0.78042, 0.672385, 0.466351])
y_der2 = y[2:] - 2 * y[1:-1] + y[0:-2]

As you can see, second derivative is negative almost everywhere, except for 1-2 points in-between, where the ball bounces and experiences big positive acceleration. So each section where it is negative corresponds to free fall. Throw away one point at each side of such section for safety.

For free fall, Y is a quadratic function in time (or in X, which practically should be the same).

Y = -0.5 * g * t^2 + b * t + c

You can then apply polynomial regression to find the quadratic equation for Y.

Then, for each parabolic-motion section, find the theoretically predicted speed at bounce point. For each bounce point, the ratio between speeds in adjacent parabolic sections is the restitution coefficient. You can verify that it is approximately constant (maybe, for debugging, also verify another theoretically predicted equality: approximately equal X-points of bounce for adjacent parabolic sections).

This method fill fail when the ball bounces too fast, so that there are not enough measurements between the bounce points to apply polynomial regression. You should stop the calculations when such a condition appears.

P.S. I think your motion model should also have vx = e * bx when the ball bounces. Not sure how important it is for prediction accuracy. Anyway, the method I described doesn't depend on this detail.

@anatolyg's method is the correct first step. After that,

• apply find_peaks;
• perform a first piecewise polynomial regression: first-degree in x and second-degree in y; then
• as a further step that I don't demonstrate, you need to enforce that the boundary positions match between each segment and then do an end-to-end fit where the physical parameters are shared across the whole dataset.

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal


def fit(
    t: np.ndarray,
    x: np.ndarray,
    y: np.ndarray,
) -> tuple[
    list[np.polynomial.Polynomial],  # x(t)
    list[np.polynomial.Polynomial],  # y(t)
]:
    # Second-order differential. This assumes a monotonic t.
    d2ydt2 = np.gradient(y, 2)

    # Boundary conditions for each bounce segment
    bounce_idx, props = scipy.signal.find_peaks(d2ydt2)
    left_indices = (0, *bounce_idx)
    right_indices = (*bounce_idx, len(t))

    # Boolean arrays selecting each segment
    segment_predicates = [
        (d2ydt2 < 0)   # Must be falling
        & (t >= left)  # Must be within second-order peaks
        & (t < right)
        for left, right in zip(left_indices, right_indices)
    ]
    x_polys = [
        np.polynomial.polynomial.Polynomial.fit(
            x=t[predicate], y=x[predicate], symbol='t', deg=1,
        )
        for predicate in segment_predicates
    ]
    y_polys = [
        np.polynomial.polynomial.Polynomial.fit(
            x=t[predicate], y=y[predicate], symbol='t', deg=2,
        )
        for predicate in segment_predicates
    ]
    return x_polys, y_polys


def dump(
    x_polys: list[np.polynomial.Polynomial],
    y_polys: list[np.polynomial.Polynomial],
) -> None:
    for i, (xp, yp) in enumerate(zip(x_polys, y_polys)):
        print(f'Bounce {i}:')
        print(f'  x={xp}')
        print(f'  y={yp}')


def plot(
    t: np.ndarray,
    x: np.ndarray,
    y: np.ndarray,
    x_polys: list[np.polynomial.Polynomial],
    y_polys: list[np.polynomial.Polynomial],
) -> plt.Figure:
    fig, ax = plt.subplots()

    ax.scatter(x, y, marker='+', label='experiment')

    tfine = np.linspace(start=t[0], stop=t[-1], num=201)
    for xp, yp in zip(x_polys, y_polys):
        xt = xp(tfine)
        yt = yp(tfine)
        use = yt >= 0
        ax.plot(xt[use], yt[use])

    return fig


def main() -> None:
    x = np.array((
        7.410000e-03, 9.591000e-02, 2.844100e-01, 5.729100e-01, 9.614100e-01,
        1.449910e+00, 2.038410e+00, 2.726910e+00, 3.373700e+00, 4.040770e+00,
        4.800040e+00, 5.577610e+00, 6.355180e+00, 7.132750e+00, 7.910320e+00,
        8.687900e+00, 9.465470e+00, 1.020976e+01, 1.092333e+01, 1.163690e+01,
        1.235047e+01, 1.306404e+01, 1.377762e+01, 1.449119e+01,
    ))
    y = np.array((
        2.991964, 2.903274, 2.716584, 2.431894, 2.049204, 1.568514, 0.989824, 0.313134,
        0.311512, 0.646741, 0.88397 , 1.0232  , 1.064429, 1.007658, 0.852887, 0.600116,
        0.249345, 0.232557, 0.516523, 0.702488, 0.790454, 0.78042 , 0.672385, 0.466351,
    ))
    t = np.arange(len(x), dtype=x.dtype)
    x_polys, y_polys = fit(t=t, x=x, y=y)
    dump(x_polys=x_polys, y_polys=y_polys)
    plot(t=t, x=x, y=y, x_polys=x_polys, y_polys=y_polys)
    plt.show()


if __name__ == '__main__':
    main()

Bounce 0:
  x=1.01716 + 1.35975·(-1.0 + 0.28571429t)
  y=2.252799 - 1.339415·(-1.0 + 0.28571429t) - 0.60025·(-1.0 + 0.28571429t)²
Bounce 1:
  x=7.910324 + 1.555146·(-7.0 + 0.5t)
  y=0.852887 - 0.407542·(-7.0 + 0.5t) - 0.196·(-7.0 + 0.5t)²
Bounce 2:
  x=13.77761667 + 0.713575·(-22.0 + t)
  y=0.672385 - 0.1570345·(-22.0 + t) - 0.0489995·(-22.0 + t)²