AI for Mechanical Engineering: Dynamics
By Keonhyeok Park
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.preprocessing.sequence import pad_sequences
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Input, RepeatVector
import os
1. Projectile Motion (Trajectory Prediction)¶
Explanation of Equations Used¶
When solving the projectile motion with air resistance, several key equations and principles are used to model the motion accurately. Below are the detailed descriptions of these equations:
Decomposing Initial Velocity¶
The initial velocity $v_0$ is decomposed into horizontal ($ v_{x0} $) and vertical ($ v_{y0} $) components based on the launch angle $ \theta $:
$$v_{x0} = v_0 cos(\theta)$$
$$v_{y0} = v_0 sin(\theta)$$
Air Resistance (Drag Force)¶
The drag force $ F_d $ acting on the projectile is given by Newton's drag law, which states that the drag force is proportional to the square of the velocity. The drag force can be calculated as:
$$F_d = \frac{1}{2} C_d \rho A v^2$$
where:
- $ C_d $ is the drag coefficient
- $ \rho $ is the air density
- $ A $ is the cross-sectional area of the projectile
- $ v $ is the instantaneous velocity of the projectile
Acceleration Components¶
The horizontal ($ a_x $) and vertical ($ a_y $) components of acceleration due to drag and gravity are given by:
$$\begin{align*} a_x &= -\frac{F_d}{m} \cdot \frac{v_x}{v} \\ a_y &= -g - \frac{F_d}{m} \cdot \frac{v_y}{v} \end{align*}$$
where:
- $ g $ is the acceleration due to gravity (9.81 m/s²)
- $ m $ is the mass of the projectile
- $ v_x $ and $ v_y $ are the horizontal and vertical components of the velocity, respectively
Solving Differential Equations¶
The projectile motion is governed by the second-order differential equations derived from Newton's second law. These equations are solved numerically using the odeint function from the scipy library. The state vector $[x, y, v_x, v_y]$ is updated at each time step.
Stopping Condition¶
The simulation time is defined as a range from 0 to a maximum time, divided into a number of time steps (e.g., 10,000). The simulation stops when the vertical position ($ y $) becomes less than or equal to zero, indicating the projectile has hit the ground. The index where this occurs is found and the data up to this point is retained:
$$ \text{ground_hit_index} = \text{np.where}(y < 0)[0][0] $$
By using these equations and principles, we can model the trajectory of a projectile under the influence of gravity and air resistance, providing a realistic simulation of its motion.
In [ ]:
def projectile_derivatives(state, t, mass, drag_coefficient, area, air_density):
x, y, vx, vy = state
# Calculate speed
v = np.sqrt(vx**2 + vy**2)
# Calculate drag force
Fd = 0.5 * drag_coefficient * air_density * area * v**2
# Calculate acceleration
ax = -Fd * (vx / v) / mass
ay = -9.81 - (Fd * (vy / v) / mass)
return [vx, vy, ax, ay]
def projectile_trajectory(v0, theta, mass, drag_coefficient, area, air_density, num_points=1000):
theta_rad = np.radians(theta) # Convert angle to radians
vx0 = v0 * np.cos(theta_rad) # Initial velocity in x-direction
vy0 = v0 * np.sin(theta_rad) # Initial velocity in y-direction
# Initial state: [x, y, vx, vy]
initial_state = [0, 0, vx0, vy0]
# Time parameters
dt = 0.01 # Time step in seconds
t = 0 # Initial time
# Lists to store trajectory points
x_points = []
y_points = []
# Initial state
state = initial_state
# Simulation loop
while state[1] >= 0:
# Append current position to trajectory
x_points.append(state[0])
y_points.append(state[1])
# Compute derivatives
derivatives = projectile_derivatives(state, t, mass, drag_coefficient, area, air_density)
# Update state using Euler's method
state = [state[i] + dt * derivatives[i] for i in range(4)]
# Update time
t += dt
# Convert lists to numpy arrays
x_points = np.array(x_points)
y_points = np.array(y_points)
# Interpolate to get fixed number of points
distance = np.linspace(0, x_points[-1], num_points)
y_interpolated = np.interp(distance, x_points, y_points)
return distance, y_interpolated
v0 = 1000 # Initial speed in m/s
theta = 80 # Launch angle in degrees
mass = 0.145 # Mass of the projectile in kg (e.g., a baseball)
drag_coefficient = 0.47 # Drag coefficient (typical for a sphere)
area = 0.0042 # Cross-sectional area in m^2 (typical for a baseball)
air_density = 1.225 # Air density in kg/m^3 (at sea level)
x, y = projectile_trajectory(v0, theta, mass, drag_coefficient, area, air_density)
# Plotting the trajectory
plt.figure(figsize=(10, 10))
plt.subplot(211)
plt.plot(x, y, '.')
plt.title("Projectile Trajectory with Air Resistance (sampling)")
plt.xlabel("Distance (m)")
plt.ylabel("Height (m)")
plt.grid(True)
plt.subplot(212)
plt.plot(x, y)
plt.title("Projectile Trajectory with Air Resistance")
plt.xlabel("Distance (m)")
plt.ylabel("Height (m)")
plt.grid(True)
plt.show()
Lab 1: Projectile Motion with Air Resistance¶
- Input: Initial angle ($\theta$), initial velocity ($v$)
- Output: Trajectory
Projectile Data Generation¶
- 1,000 dataset
- Angles: 10 to 80 degrees
- Velocities: 50 to 1000 m/s
In [ ]:
def generate_projectile_data(n_samples=1000):
angles = np.linspace(10, 80, n_samples)
velocities = np.linspace(50, 1000, n_samples)
data = []
for v0, theta in zip(velocities, angles):
x, y = projectile_trajectory(v0, theta, mass, drag_coefficient, area, air_density)
trajectory = np.vstack((x, y)).T
data.append((v0, theta, trajectory))
return data
data = generate_projectile_data()
In [ ]:
random_idx = np.array([0,100,500])
plt.figure(figsize = (15,8))
plt.title('Velocity: {:.2f} / {:.2f} / {:.2f}\nAngle: {:.2f} / {:.2f} / {:.2f}'.format(data[random_idx[0]][0],data[random_idx[1]][0],data[random_idx[2]][0],
data[random_idx[0]][1],data[random_idx[1]][1],data[random_idx[2]][1]),
fontsize=20)
for idx in random_idx:
plt.plot(data[idx][2][:,0],data[idx][2][:,1], label = idx, linewidth=3, linestyle='-')
plt.legend(title='Random Index',title_fontsize = 15, fontsize = 15)
Out[Â ]:
- Preprocess the dataset
In [ ]:
def preprocess_data(data):
max_len = max(len(traj[2]) for traj in data)
X = []
y = []
for v0, theta, trajectory in data:
input_features = np.array([v0, theta])
trajectory = np.array(trajectory)
X.append(input_features) # Repeat input features to match trajectory length
y.append(trajectory)
X = np.array(X)
y = np.array(y)
input_scaler = MinMaxScaler(feature_range=(0, 1))
output_scaler = MinMaxScaler(feature_range=(0, 1))
# Flatten the input and output arrays for scaling
X_flattened = X.reshape(-1, X.shape[-1])
y_flattened = y.reshape(-1, y.shape[-1])
# Fit and transform the scalers
X_scaled_flattened = input_scaler.fit_transform(X_flattened)
y_scaled_flattened = output_scaler.fit_transform(y_flattened)
# Reshape back to the original shape
X_scaled = X_scaled_flattened.reshape(X.shape)
y_scaled = y_scaled_flattened.reshape(y.shape)
return X_scaled, y_scaled, input_scaler, output_scaler
X_scaled, y_scaled, input_scaler, output_scaler = preprocess_data(data)
# Display shapes to verify
print("X_scaled shape:", X_scaled.shape)
print("y_scaled shape:", y_scaled.shape)
LSTM Model Architecture¶
In [ ]:
model = Sequential()
model.add(Dense(128, activation='relu', input_shape=(2,))) # (batch_size, 2) -> (batch_size, 128)
model.add(RepeatVector(y_scaled.shape[1])) # (batch_size, 1000, 128)
model.add(LSTM(128, return_sequences=True)) # (batch_size, 1000, 128)
model.add(Dense(2)) # (batch_size, 1000, 2) == (batch_size, timesteps, features)
# Model compile
model.compile(optimizer='adam', loss='mse')
model.summary()
In [ ]:
history = model.fit(X_scaled, y_scaled, epochs=50, batch_size=32, validation_split=0.2)
In [ ]:
plt.plot(history.history['loss'], label='train loss')
plt.plot(history.history['val_loss'], label='val loss')
plt.legend()
plt.show()
In [ ]:
v0_test = 300
theta_test = 45
input_features = np.array([v0_test, theta_test]).reshape(1, -1)
input_scaled = input_scaler.transform(input_features)
prediction_scaled = model.predict(input_scaled)
predicted_trajectory = output_scaler.inverse_transform(prediction_scaled.reshape(-1, 2))
true_distance, true_height = projectile_trajectory(v0_test, theta_test, mass, drag_coefficient, area, air_density)
plt.plot(predicted_trajectory[:, 0], predicted_trajectory[:, 1], label='Predicted Trajectory')
plt.plot(true_distance, true_height, label='True Trajectory', linestyle='dashed')
plt.xlabel('Distance (m)')
plt.ylabel('Height (m)')
plt.title('Predicted vs True Projectile Trajectory')
plt.legend()
plt.show()
Ineffective Cases¶
- Input: Initial angles, initial velocities
- Output: Trajectory
[Reasons for Ineffectiveness]
- Predicting multiple future time points based only on the current state may lead to insufficient information, especially when the model needs to learn nonlinearities and complex interactions.
- Long-term predictions based solely on initial conditions are very challenging. Small prediction errors can accumulate over time, leading to significant inaccuracies.
Lab 2: Projectile Motion with Air Resistance¶
- Input: Partial Trajectory (10 timesteps)
- Output: Next Position (1 timestep)
Data Generation¶
In [ ]:
def projectile_trajectory(v0, theta, mass, drag_coefficient, area, air_density, num_points=100):
theta_rad = np.radians(theta) # Convert angle to radians
vx0 = v0 * np.cos(theta_rad) # Initial velocity in x-direction
vy0 = v0 * np.sin(theta_rad) # Initial velocity in y-direction
# Initial state: [x, y, vx, vy]
initial_state = [0, 0, vx0, vy0]
# Time parameters
dt = 0.01 # Time step in seconds
t = 0 # Initial time
# Lists to store trajectory points
x_points = []
y_points = []
# Initial state
state = initial_state
# Simulation loop
while state[1] >= 0:
# Append current position to trajectory
x_points.append(state[0])
y_points.append(state[1])
# Compute derivatives
derivatives = projectile_derivatives(state, t, mass, drag_coefficient, area, air_density)
# Update state using Euler's method
state = [state[i] + dt * derivatives[i] for i in range(4)]
# Update time
t += dt
# Convert lists to numpy arrays
x_points = np.array(x_points)
y_points = np.array(y_points)
# Interpolate to get fixed number of points
distance = np.linspace(0, x_points[-1], num_points)
y_interpolated = np.interp(distance, x_points, y_points)
return distance, y_interpolated
In [ ]:
train_window = 10
def create_inout_sequences(input_data, tw):
inout_seq = []
L = len(input_data)
for i in range(90):
train_seq = input_data[i:i+tw]
train_label = input_data[i+tw:i+tw+1]
inout_seq.append((train_seq ,train_label))
return inout_seq
def make_inout_seq(angle, speed):
X, Y = projectile_trajectory(speed, angle, mass, drag_coefficient, area, air_density)
x = np.zeros((len(X), 2))
for i in range(len(X)):
x[i, 0] = X[i]
x[i, 1] = Y[i]
scaler = MinMaxScaler(feature_range=(0, 1))
train_data_normalized = scaler.fit_transform(x.reshape(-1, 2))
train_window = 10
train_inout_seq = create_inout_sequences(train_data_normalized.reshape(len(X), 2), train_window)
return scaler, train_inout_seq
scaler, train_inout_seq = make_inout_seq(15, 500)
In [ ]:
X_train = np.array([seq[0] for seq in train_inout_seq])
y_train = np.array([seq[1] for seq in train_inout_seq]).squeeze()
print(X_train.shape, y_train.shape)
In [ ]:
plt.figure(figsize=(12,5))
for i in range(15):
scaler_test, test_inout_seq = make_inout_seq(5*i, 500)
X_test = np.array([seq[0] for seq in test_inout_seq])
y_test = np.array([seq[1] for seq in test_inout_seq])
X_test_traj = scaler_test.inverse_transform(X_test[0].reshape(-1,2))
y_test_traj = scaler_test.inverse_transform(y_test.reshape(-1,2))
plt.plot(X_test_traj[:,0], X_test_traj[:,1])
plt.plot(y_test_traj[:,0], y_test_traj[:,1], label = 'actual', color = 'blue')
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(), fontsize = 15)
plt.xlim([0,400])
plt.ylim([0,300])
plt.title("Projectile Trajectory with Air Resistance", fontsize = 15)
plt.xlabel("Distance (m)", fontsize = 15)
plt.ylabel("Height (m)", fontsize = 15)
Out[Â ]:
LSTM Model Architecture¶
In [ ]:
# Build LSTM model with Keras
model = Sequential()
model.add(LSTM(100, input_shape=(train_window, 2), return_sequences=True))
model.add(LSTM(100))
model.add(Dense(100))
model.add(Dense(2))
model.compile(optimizer='adam', loss='mse')
model.summary()
In [ ]:
history = model.fit(X_train, y_train, epochs=400, batch_size=32)
In [ ]:
plt.plot(history.history['loss'], label='train loss')
# plt.plot(history.history['val_loss'], label='val loss')
plt.legend()
plt.show()
Rollout Prediction¶
In [ ]:
plt.figure(figsize=(12,5))
for i in range(15):
scaler_test, test_inout_seq = make_inout_seq(5*i, 500)
X_test = np.array([seq[0] for seq in test_inout_seq])
y_test = np.array([seq[1] for seq in test_inout_seq])
test_seq = X_test[0]
for i in range(90):
pred = model.predict(test_seq[i:][np.newaxis], verbose = 0)
test_seq = np.concatenate([test_seq, pred], axis = 0).squeeze()
pred_traj = scaler_test.inverse_transform(test_seq.reshape(-1,2))
pred_traj = pred_traj[pred_traj[:,1]>=0]
X_test_traj = scaler_test.inverse_transform(X_test[0].reshape(-1,2))
y_test_traj = scaler_test.inverse_transform(y_test.reshape(-1,2))
plt.plot(X_test_traj[:,0], X_test_traj[:,1])
plt.plot(y_test_traj[:,0], y_test_traj[:,1], label = 'actual', color = 'blue')
plt.plot(np.concatenate([X_test_traj[train_window-1,0][np.newaxis], pred_traj[:,0]], axis = 0),
np.concatenate([X_test_traj[train_window-1,1][np.newaxis], pred_traj[:,1]], axis = 0), label = 'pred', color = 'red')
# plt.vlines(X_test_traj[:train_window,0][-1],0,np.max(y_test_traj[:,1])*1.1, color = 'red', linestyle='--')
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(), fontsize = 15)
plt.title("Projectile Trajectory with Air Resistance", fontsize = 15)
plt.xlabel("Distance (m)", fontsize = 15)
plt.ylabel("Height (m)", fontsize = 15)
plt.grid(True)
plt.xlim([0,400])
plt.ylim([0,300])
Out[Â ]: