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[ ]:
Lab 3: Projectile Motion with Air Resistance (PINN)¶
- Physics-informed Neural Network (PINN)
In [ ]:
g = 9.81
C_d = 0.47
rho = 1.225
A = 0.0042
m = 0.145
class PINN(tf.keras.Model):
def __init__(self):
super(PINN, self).__init__()
self.dense1 = tf.keras.layers.Dense(50, activation='tanh')
self.dense2 = tf.keras.layers.Dense(50, activation='tanh')
self.dense3 = tf.keras.layers.Dense(50, activation='tanh')
self.dense4 = tf.keras.layers.Dense(2) # Output: [x, y]
def call(self, t):
x = self.dense1(t)
x = self.dense2(x)
x = self.dense3(x)
return self.dense4(x)
def loss(model, t, x0, y0, vx0, vy0):
with tf.GradientTape(persistent=True) as tape:
tape.watch(t)
pred = model(t)*300
x, y = pred[:, 0:1], pred[:, 1:2]
x_t = tape.gradient(x, t)
y_t = tape.gradient(y, t)
x_tt = tape.gradient(x_t, t)
y_tt = tape.gradient(y_t, t)
# PDE residuals
v = tf.sqrt(x_t**2 + y_t**2)
f_x = x_tt + (0.5 * C_d * rho * A * x_t * v) / m
f_y = y_tt + g + (0.5 * C_d * rho * A * y_t * v) / m
ic_x = x[0] - x0
ic_y = y[0] - y0
ic_vx = x_t[0] - vx0
ic_vy = y_t[0] - vy0
return tf.reduce_mean(tf.square(f_x)) + \
tf.reduce_mean(tf.square(f_y)) + \
tf.square(ic_x) + \
tf.square(ic_y) + \
tf.square(ic_vx) + \
tf.square(ic_vy)
v0 = 500.0
theta = np.pi / 12 * 3
x0 = 0.0
y0 = 0.0
vx0 = v0 * np.cos(theta)
vy0 = v0 * np.sin(theta)
t = np.linspace(0, 15, 1000).reshape(-1, 1)
t = tf.convert_to_tensor(t, dtype=tf.float32)
In [ ]:
model.summary()
In [ ]:
model = PINN()
optimizer = tf.keras.optimizers.Adam(learning_rate=0.001)
epochs = 10000
for epoch in range(epochs):
with tf.GradientTape() as tape:
current_loss = loss(model, t, x0, y0, vx0, vy0)
gradients = tape.gradient(current_loss, model.trainable_variables)
optimizer.apply_gradients(zip(gradients, model.trainable_variables))
if epoch % 500 == 0:
print(f"Epoch {epoch}: Loss = {current_loss.numpy()}")
Evaluation¶
In [ ]:
t_test = np.linspace(0, 20, 1000).reshape(-1, 1)
t_test = tf.convert_to_tensor(t_test, dtype=tf.float32)
pred = model(t_test)*300
x_pred, y_pred = pred[:, 0].numpy(), pred[:, 1].numpy()
valid_indices = y_pred >= 0
x_pred = x_pred[valid_indices]
y_pred = y_pred[valid_indices]
def true_trajectory_with_drag(dt=0.01):
t = 0
x, y = [0], [0]
vx, vy = vx0, vy0
while y[-1] >= 0:
v = np.sqrt(vx**2 + vy**2)
ax = - (0.5 * C_d * rho * A * v * vx) / m
ay = - g - (0.5 * C_d * rho * A * v * vy) / m
vx += ax * dt
vy += ay * dt
x.append(x[-1] + vx * dt)
y.append(y[-1] + vy * dt)
t += dt
print(t)
return np.array(x), np.array(y)
x_true, y_true = true_trajectory_with_drag()
plt.figure(figsize=(12,5))
plt.plot(x_true, y_true, 'blue', label='Actual')
plt.plot(x_pred, y_pred, label='Predicted', color='red')
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 (PINN)", fontsize = 15)
plt.xlabel("Distance (m)", fontsize = 15)
plt.ylabel("Height (m)", fontsize = 15)
plt.grid(True)
plt.xlim([0,400])
plt.ylim([0,310])
Out[ ]:
Spring Mass System¶
Lab 1: Coupled Spring Mass ODE System¶
The provided diagram depicts a system consisting of three masses ($m_1$, $m_2$, and $m_3$) connected by springs in a linear arrangement. The springs are connected to rigid walls at both ends. To analyze the dynamics of this system, we need to derive the Partial Differential Equations (PDEs) governing the motion of each mass.
- $m_1$, $m_2$, and $m_3$ are masses connected in series with springs.
- The springs have stiffness constants $k$.
- The system is attached to rigid walls on both ends.
- The displacements of the masses are denoted as $x_1(t)$, $x_2(t)$, and $x_3(t)$ respectively.
Equations of Motion
$$ \begin{cases} m_1 \ddot{x}_1(t) &= -k_1 x_1(t) + k_2 \left( x_2(t) - x_1(t) \right), \\ m_2 \ddot{x}_2(t) &= -k_2 \left( x_2(t) - x_1(t) \right) + k_3 \left( x_3(t) - x_2(t) \right), \\ m_3 \ddot{x}_3(t) &= -k_3 \left( x_3(t) - x_2(t) \right) - k_4 x_3(t). \end{cases} $$
Parameters
$$ \begin{align*} [m_1, m_2, m_3] &= [1.0, 1.0, 1.0], \\ [k_1, k_2, k_3, k_4] &= [2.0, 1.0, 1.0, 2.0], \end{align*} $$
Initial Conditions
$$ \begin{align*} [x_1(0), x_2(0), x_3(0)] &= [1.0, 0.5, 0.0], \\ [\dot{x}_1(0), \dot{x}_2(0), \dot{x}_3(0)] &= [0.0, 0.0, 0.0]. \end{align*} $$
- Input: Partial Trajectory (10 timesteps)
- Output: Next Position (1 timestep)
- Train (0 seconds to 10 seconds)
- Test (0 seconds to 20 seconds)
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Reshape, Concatenate
import matplotlib.animation as animation
from sklearn.preprocessing import MinMaxScaler
In [ ]:
masses = np.array([1.0, 1.0, 1.0])
spring_constants = np.array([2.0, 1.0, 1.0, 2.0])
n_masses = len(masses)
n_springs = len(spring_constants)
In [ ]:
initial_positions = np.array([1,0.5,0])
initial_velocities = np.zeros(n_masses)
initial_conditions = np.concatenate([initial_positions, initial_velocities])
t_span = [0, 10] # 시간 범위
t_eval = np.linspace(t_span[0], t_span[1], 200)
print(initial_conditions, t_eval.shape)
In [ ]:
t_span_test = [0, 20]
t_eval_test = np.linspace(t_span_test[0], t_span_test[1], 400)
In [ ]:
# motion equation
def motion(t, y):
positions = y[:n_masses]
velocities = y[n_masses:]
accelerations = np.zeros(n_masses)
for i in range(n_masses):
if i > 0:
accelerations[i] -= spring_constants[i] * (positions[i] - positions[i-1]) / masses[i]
if i < n_masses - 1:
accelerations[i] += spring_constants[i+1] * (positions[i+1] - positions[i]) / masses[i]
accelerations[0] -= spring_constants[0] * positions[0] / masses[0]
accelerations[-1] += spring_constants[-1] * (0 - positions[-1]) / masses[-1]
return np.concatenate([velocities, accelerations])
In [ ]:
solution = solve_ivp(motion, t_span, initial_conditions, t_eval=t_eval, method='RK45')
solution_test = solve_ivp(motion, t_span_test, initial_conditions, t_eval=t_eval_test, method='RK45')
In [ ]:
fig, axs = plt.subplots(1, 1, figsize=(10, 5))
# Position graphs for each mass
for i in range(n_masses):
axs.plot(t_eval, solution.y[i], label=f'Position of Mass {i+1}')
axs.set_xlabel('Time (s)', fontsize = 15)
axs.set_ylabel('Position', fontsize = 15)
axs.set_title('Position vs Time for Each Mass', fontsize = 15)
axs.vlines(10, -1, 1.1, color='red',linestyle='--')
axs.legend(fontsize = 15, loc='upper right')
axs.set_xlim([0,20])
plt.xlim([0,20])
plt.tight_layout()
plt.show()
In [ ]:
############################ TEST data ################################
fig, axs = plt.subplots(1, 1, figsize=(10, 5))
# Position graphs for each mass
for i in range(n_masses):
axs.plot(t_eval_test, solution_test.y[i], label=f'Position of Mass {i+1}')
axs.set_xlabel('Time (s)', fontsize = 15)
axs.set_ylabel('Position', fontsize = 15)
axs.set_title('Position vs Time for Each Mass', fontsize = 15)
axs.legend(fontsize = 15, loc='upper right')
axs.vlines(10, -1, 1.1, color='red',linestyle='--')
plt.tight_layout()
plt.show()
Data Generation¶
In [ ]:
def create_inout_sequences(input_data, tw):
inout_seq = []
L = len(input_data)
for i in range(L-tw):
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
data = solution.y[:3].T
scaler = MinMaxScaler(feature_range=(-1,1))
train_data_normalized = scaler.fit_transform(data)
train_window = 10
train_inout_seq = create_inout_sequences(train_data_normalized, train_window)
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 [ ]:
fig, axs = plt.subplots(1, 1, figsize=(10, 4))
# Position graphs for each mass
for i in range(n_masses):
axs.plot(t_eval, solution.y[i], alpha = 0.3, color = 'black')
axs.plot(t_eval[140:150], solution.y[i][140:150], label=f'Position of Mass {i+1}')
axs.set_xlabel('Time (s)', fontsize = 15)
axs.set_ylabel('Position', fontsize = 15)
axs.set_title('Position vs Time for Each Mass', fontsize = 15)
axs.legend(fontsize = 15, loc='upper right')
axs.vlines(7, -1, 1.1, color='red',linestyle='--')
axs.vlines(7.5, -1, 1.1, color='red',linestyle='--')
# # Velocity graphs for each mass
# for i in range(n_masses):
# axs[1].plot(t_eval, solution.y[i + n_masses], label=f'Velocity of Mass {i+1}')
# axs[1].set_xlabel('Time (s)', fontsize = 15)
# axs[1].set_ylabel('Velocity', fontsize = 15)
# axs[1].set_title('Velocity vs Time for Each Mass', fontsize = 15)
# axs[1].legend(fontsize = 15, loc='upper right')
plt.xlim([0,11])
plt.ylim([-1,1.1])
plt.tight_layout()
plt.show()
LSTM Model Architecture¶
In [ ]:
# Build LSTM model with Keras
model = Sequential()
model.add(LSTM(100, input_shape=(None, 3), return_sequences=True))
model.add(LSTM(100, activation = 'tanh'))
model.add(Dense(100, activation = 'tanh'))
model.add(Dense(3, activation = 'tanh'))
model.compile(optimizer='adam', loss='mse')
model.summary()
In [ ]:
history = model.fit(X_train, y_train, epochs=10, 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()
In [ ]:
data_test = solution_test.y[:3].T
test_data_normalized = scaler.fit_transform(data_test)
test_seq = test_data_normalized[:200]
for i in range(200):
pred = model.predict(test_seq[i:][np.newaxis], verbose = 0)
test_seq = np.concatenate([test_seq, pred], axis = 0).squeeze()
pred_traj = scaler.inverse_transform(test_seq.reshape(-1,3))
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(10,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:,0], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:,0],
label='Predicted', linestyle='--', linewidth =3, color = 'red')
plt.plot(t_eval, solution.y[0], linewidth = 3)
plt.title('Actual vs Predicted Positions (Mass 1)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(10,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:,1], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:,1],
label='Predicted', linestyle='--', linewidth =3, color = 'red')
plt.plot(t_eval, solution.y[1], linewidth = 3, color='tab:orange')
plt.title('Actual vs Predicted Positions (Mass 2)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(10,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:,2], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:,2],
label='Predicted', linestyle='--', linewidth =3, color = 'red')
plt.plot(t_eval, solution.y[2], linewidth = 3, color='tab:green')
plt.title('Actual vs Predicted Positions (Mass 3)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
Lab 2: Coupled Spring Mass ODE System (PINN)¶
In [ ]:
# !pip install deepxde/
In [ ]:
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
from deepxde.backend.set_default_backend import set_default_backend
set_default_backend("pytorch")
# Masses
m1 = 1.0
m2 = 1.0
m3 = 1.0
# Spring constants
k1 = 2.0
k2 = 1.0
k3 = 1.0
k4 = 2.0
# Initial conditions
# Initial displacements (from equilibrium positions)
d1_0 = 1
d2_0 = 0
d3_0 = 0
# Initial velocities
v1_0 = 0.0
v2_0 = 0.0
v3_0 = 0.0
# Maximum time to simulate
t_max = 5
def dy(t, x):
return dde.grad.jacobian(x, t)
def pde(t, x):
# mass 1 location
x_1 = x[:, 0:1] # Use slicing to extract the column as a tensor
# mass 2 location
x_2 = x[:, 1:2] # Use slicing to extract the column as a tensor
# mass 3 location
x_3 = x[:, 2:3] # Use slicing to extract the column as a tensor
dx1_tt = dde.grad.hessian(x, t, i=0,j=0,component=0)
dx2_tt = dde.grad.hessian(x, t, i=0,j=0,component=1)
dx3_tt = dde.grad.hessian(x, t, i=0,j=0,component=2)
pde1 = m1 * dx1_tt + k1 * (x_1) - k2 * (x_2 - x_1)
pde2 = m2 * dx2_tt + k2 * (x_2 - x_1) - k3 * (x_3 - x_2)
pde3 = m3 * dx3_tt + k3 * (x_3 - x_2) - k4 * (x_3)
return [pde1, pde2, pde3]
def boundary_init(t, on_boundary):
return on_boundary and np.isclose(t[0], 0)
geom = dde.geometry.Interval(0, t_max)
init_d1 = dde.icbc.PointSetBC(np.array([0]), np.array([d1_0]).reshape(-1, 1), component=0)
init_d2 = dde.icbc.PointSetBC(np.array([0]), np.array([d2_0]).reshape(-1, 1), component=1)
init_d3 = dde.icbc.PointSetBC(np.array([0]), np.array([d3_0]).reshape(-1, 1), component=2)
init_v1 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y[:, 0:1]), boundary_init)
init_v2 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y[:, 1:2]), boundary_init)
init_v3 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y[:, 2:3]), boundary_init)
data = dde.data.PDE(
geom,
pde,
[init_d1, init_d2, init_d3, init_v1, init_v2, init_v3],
num_domain=2000,
num_boundary=100,
num_test=1000
)
layer_size = [1] + [20] * 3 + [3]
activation = "tanh"
initializer = "Glorot uniform"
In [ ]:
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr=1e-3)
In [ ]:
losshistory, train_state = model.train(epochs=5000)
dde.saveplot(losshistory, train_state, issave=False, isplot=False)
In [ ]:
dde.optimizers.config.set_LBFGS_options(maxiter=500)
model.compile("L-BFGS")
losshistory, train_state = model.train()
In [ ]:
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
In [ ]:
from scipy.integrate import odeint
# Define the differential equations for the coupled spring-mass system for ODE solver
def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.
Arguments:
w : vector of the state variables:
w = [x1,y1,x2,y2,x3,y3]
t : time
p : vector of the parameters:
p = [m1,m2,m3,k1,k2,k3,k4,L1,L2,L3,L4]
"""
x1, y1, x2, y2, x3, y3 = w
m1, m2, m3, k1, k2, k3, k4 = p
# Create f = (x1',y1',x2',y2',x3',y3'):
f = [y1,
(-k1 * (x1) + k2 * (x2 - x1)) / m1,
y2,
(-k2 * (x2 - x1) + k3 * (x3 - x2)) / m2,
y3,
(-k3 * (x3 - x2) - k4 * (x3)) / m3]
return f
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = t_max
numpoints = 250
t = geom.random_points(5000)
t[:,0].sort()
# Pack up the parameters and initial conditions:
p = [m1, m2, m3, k1, k2, k3, k4]
w0 = [d1_0, v1_0, d2_0, v2_0, d3_0, v3_0]
# Call the ODE solver.
wsol = odeint(vectorfield, w0, t[:,0], args=(p,), atol=abserr, rtol=relerr)
In [ ]:
# Predict using the trained model
result = model.predict(t)
usol1 = np.array(result[:, 0])
usol2 = np.array(result[:, 1])
usol3 = np.array(result[:, 2])
# Plot the results
lw = 2
plt.figure(figsize=(12, 6))
plt.plot(t, wsol[:, 0], alpha=1, label='Mass1 GT', c='r', lw=lw)
plt.plot(t, wsol[:, 2], alpha=1, label='Mass2 GT', c='b', lw=lw)
plt.plot(t, wsol[:, 4], alpha=1, label='Mass3 GT', c='g', lw=lw)
plt.plot(t, usol1, alpha=1, label='Mass1 Prediction', linestyle='dashed', c='k', lw=lw)
plt.plot(t, usol2, alpha=1, label='Mass2 Prediction', linestyle='dashed', c='m', lw=lw)
plt.plot(t, usol3, alpha=1, label='Mass3 Prediction', linestyle='dashed', c='y', lw=lw)
plt.legend(fontsize=15)
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)
plt.xlabel('Time (s)', fontsize=30)
plt.ylabel('Displacement (m)', fontsize=30)
plt.show()
Lab 3: LSTM with Physics¶
- Input: Partial Trajectory (10 timesteps)
- Output: Next Positions (3 timestep)
- Train (0 seconds to 5 seconds)
- Test (0 seconds to 10 seconds)
In [ ]:
masses = np.array([1.0, 1.0, 1.0])
spring_constants = np.array([2.0, 1.0, 1.0, 2.0])
n_masses = len(masses)
n_springs = len(spring_constants)
In [ ]:
initial_positions = np.array([1,0,0])
initial_velocities = np.zeros(n_masses)
initial_conditions = np.concatenate([initial_positions, initial_velocities])
t_span = [0, 5] # 시간 범위
t_eval = np.linspace(t_span[0], t_span[1], 100)
print(initial_conditions, t_eval.shape)
In [ ]:
t_span_test = [0, 10]
t_eval_test = np.linspace(t_span_test[0], t_span_test[1], 200)
In [ ]:
# motion equation
def motion(t, y):
positions = y[:n_masses]
velocities = y[n_masses:]
accelerations = np.zeros(n_masses)
for i in range(n_masses):
if i > 0:
accelerations[i] -= spring_constants[i] * (positions[i] - positions[i-1]) / masses[i]
if i < n_masses - 1:
accelerations[i] += spring_constants[i+1] * (positions[i+1] - positions[i]) / masses[i]
accelerations[0] -= spring_constants[0] * positions[0] / masses[0]
accelerations[-1] += spring_constants[-1] * (0 - positions[-1]) / masses[-1]
return np.concatenate([velocities, accelerations])
In [ ]:
solution = solve_ivp(motion, t_span, initial_conditions, t_eval=t_eval, method='RK45')
solution_test = solve_ivp(motion, t_span_test, initial_conditions, t_eval=t_eval_test, method='RK45')
Data Generation¶
In [ ]:
def create_inout_sequences(input_data, tw):
inout_seq = []
L = len(input_data)
for i in range(L-tw-2):
train_seq = input_data[i:i+tw]
train_label = input_data[i+tw:i+tw+3]
inout_seq.append((train_seq ,train_label))
return inout_seq
data = solution.y[:3].T
scaler = MinMaxScaler(feature_range=(-1,1))
train_data_normalized = scaler.fit_transform(data)
train_window = 10
train_inout_seq = create_inout_sequences(train_data_normalized, train_window)
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)
LSTM Model Architecture¶
In [ ]:
from tensorflow.keras import backend as K
def PhysicsGuidedLoss(y_true, y_pred):
# Loss 1
mse_loss = K.mean((y_true-y_pred)**2)
# Loss 2
delta_t = t_eval[1]
a_pred = (y_pred[:,2]-2*y_pred[:,1]+y_pred[:,0]) / (delta_t**2)
x1_pred, x2_pred, x3_pred = y_pred[:, :, 0], y_pred[:, :, 1], y_pred[:, :, 2]
physics_loss_1 = a_pred[:, np.newaxis, 0] - (-spring_constants[0] * x1_pred + spring_constants[1] * (x2_pred - x1_pred)) / masses[0]
physics_loss_2 = a_pred[:, np.newaxis, 1] - (-spring_constants[1] * (x2_pred - x1_pred) + spring_constants[2] * (x3_pred - x2_pred)) / masses[1]
physics_loss_3 = a_pred[:, np.newaxis, 2] - (-spring_constants[2] * (x3_pred - x2_pred) - spring_constants[3] * x3_pred) / masses[2]
physics_loss = tf.reduce_mean(tf.square(physics_loss_1) + tf.square(physics_loss_2) + tf.square(physics_loss_3))
# Total Loss
total_loss = mse_loss + physics_loss*1e-4
return total_loss
In [ ]:
# Build LSTM model with Keras
model = Sequential()
model.add(LSTM(100, input_shape=(10, 3), return_sequences=True))
model.add(LSTM(100, activation = 'tanh'))
model.add(Dense(100, activation = 'tanh'))
model.add(Dense(9, activation = 'tanh'))
model.add(Reshape((3,3)))
model.compile(optimizer='adam', loss=PhysicsGuidedLoss)
model.summary()
In [ ]:
history = model.fit(X_train, y_train, epochs=500, 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()
In [ ]:
data_test = solution_test.y[:3].T
test_data_normalized = scaler.fit_transform(data_test)
test_seq = test_data_normalized[:10]
for i in range(200):
pred = model.predict(test_seq[3*i:][np.newaxis], verbose = 0)
test_seq = np.concatenate([test_seq, pred.squeeze()], axis = 0).squeeze()
pred_traj = scaler.inverse_transform(test_seq.reshape(-1,3))
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(0.45,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:200,0], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:200,0],
label='Predicted', linestyle='-', linewidth =3, color = 'red')
plt.plot(t_eval[:10], solution.y[0][:10], linewidth = 3)
plt.title('Actual vs Predicted Positions (Mass 1)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(0.45,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:200,1], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:200,1],
label='Predicted', linestyle='-', linewidth =3, color = 'red')
plt.plot(t_eval[:10], solution.y[1][:10], linewidth = 3, color='tab:orange')
plt.title('Actual vs Predicted Positions (Mass 2)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
In [ ]:
plt.figure(figsize=(15,7))
plt.vlines(0.45,-1,1, color = 'red', linestyle='--')
plt.plot(t_eval_test, data_test[:200,2], label='Actual', linewidth = 3, color='black',alpha=0.3)
plt.plot(t_eval_test, pred_traj[:200,2],
label='Predicted', linestyle='-', linewidth =3, color = 'red')
plt.plot(t_eval[:10], solution.y[2][:10], linewidth = 3, color='tab:green')
plt.title('Actual vs Predicted Positions (Mass 3)', fontsize = 15)
plt.xlabel('Time', fontsize = 15)
plt.ylabel('Position', fontsize = 15)
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.show()
Trajectory of Double Pendulum¶
What is double pendulum problem?
Complexitiy of the Dobule Pendulum Problem
- Nonlinearity
- The double pendulum's behavior is governed by nonlinear equations, leading to complex, unpredictable dynamics.
- Sensitivity to Initial Conditions
- Minor variations in the initial state can result in drastically different motions.
-Lack of Periodicity
- Unlike simple pendulums, the double pendulum lacks consistent cycles
- Dependence on Multiple Parameters
- The interplay of parameters like mass, length, and gravity affects the pendulum's motion.
In [ ]:
import numpy as np
from scipy.integrate import odeint
import math
from sklearn.preprocessing import MinMaxScaler
import tensorflow as tf
from matplotlib import animation, rc
from matplotlib import pyplot as plt
from IPython.display import HTML
from sklearn.metrics import mean_absolute_error, mean_squared_error
In [ ]:
# Constants
m1 = 2.0 # Mass 1
m2 = 1.0 # Mass 2
L1 = 1.5 # Length 1
L2 = 1.0 # Length 2
g = 9.8 # Acceleration due to gravity
Configuration vector $u ~=~ [\theta_1, \dot{\theta}_1, \theta_2, \dot{\theta}_2 ]$ has two angles $\theta_1$ and $\theta_2$ and the angular velocity $\dot{\theta}_1$ and $\dot{\theta}_2$ so
$$ \dot{\theta}_1 = \frac{d \theta_1}{dt} ~~~ and ~~~~ \dot{\theta}_2 = \frac{d \theta_2}{dt} $$
For the purposes of numerical solution, the equations are written in the Lagrangian formalism.
$$ \ddot{\theta}_1 ~=~ \frac{- (m_1+m_2)g sin\theta_1 + m_2 g sin\theta_2 sin\Delta - m_2 sin\Delta (L_2 \dot\theta_2^2 + L_1\dot\theta_1^2cos\Delta)}{ (m_1+m_2)L_1 - m_2L_1cos^2\Delta} $$
$$ \ddot{\theta}_2 ~=~ \frac{(m_1+m_2)(L_1\dot\theta_1^2sin\Delta - gsin\dot\theta_2 + gsin\dot\theta_1cos\Delta) + m_2 L_2 \dot\theta_2^2 sin\Delta cos\Delta}{(m_1+m_2)L_2 - m_2L_2cos^2\Delta} $$
where $m_i, L_i$ and $\theta_i$ denote the masses, lengths, and angles from the vertical of the upper $(i = 1)$ and lower $(i = 2)$ pendulum, and $\theta_1 - \theta_2 = \Delta$.
In [ ]:
def derivs(state, t):
state_derivs = np.zeros_like(state)
theta1, theta1_dot, theta2, theta2_dot = state
delta = theta1 - theta2
c = np.cos(delta)
s = np.sin(delta)
state_derivs[0] = theta1_dot # d(theta1)/dt
state_derivs[1] = ((- (m1 + m2) * g * np.sin(theta1)
+ m2 * g * np.sin(theta2) * c
- m2 * s * (L2 * theta2_dot**2 + L1 * theta1_dot**2 * c))
/ ((m1 + m2) * L1 - m2 * L1 * c**2))
state_derivs[2] = theta2_dot # d(theta2)/dt
state_derivs[3] = (((m1 + m2) * (L1 * theta1_dot**2 * s
- g * np.sin(theta2)
+ g * np.sin(theta1) * c)
+ m2 * L2 * theta2_dot**2 * s * c)
/ ((m1 + m2) * L2 - m2 * L2 * c**2))
return state_derivs
Lab 1: LSTM for Prediction¶
In [ ]:
# Simulation parameters
tfinal = 20.0 # Final time. Simulation time = 0 to tfinal.
Nt = 1000
t = np.linspace(0, tfinal, Nt)
theta1 = -80
theta2 = -150
# Hard coded initial velocities, can change to random
initial_velo_1 = 0.0
initial_velo_2 = 0.0
# Running the simulation
state = np.radians([theta1, initial_velo_1, theta2, initial_velo_2])
sol = odeint(derivs, state, t)
In [ ]:
#Generate the training set
train_window = 50
def create_inout_sequences(input_data, tw):
x = []
y = []
L = len(input_data)
for i in range(L-tw):
seq_tmp = input_data[i:i+tw]
label_tmp = input_data[i+tw:i+tw+1]
x.append(seq_tmp)
y.append(label_tmp)
return np.array(x), np.array(y)
def preprocess_sol(sol, train_window):
scaler = MinMaxScaler(feature_range=(-1, 1)) # -1 ~ 1 normalization
train_data_normalized = scaler.fit_transform(sol)
train_data_normalized = np.array(train_data_normalized)
train_X, train_y= create_inout_sequences(train_data_normalized, train_window)
return train_X, train_y # (990, 10, 4)
In [ ]:
scaler = MinMaxScaler(feature_range=(-1, 1)) # -1 ~ 1 normalization
train_data_normalized = scaler.fit_transform(sol)
train_data_normalized = np.array(train_data_normalized)
data, label= create_inout_sequences(train_data_normalized, train_window)
train_data = np.array(data).reshape((-1,train_window,4))[:500]
test_data = np.array(data).reshape((-1,train_window,4))[500:]
train_label = np.array(label).reshape((-1,4))[:500]
test_label = np.array(label).reshape((-1,4))[500:]
Model Architecture¶
In [ ]:
n_lstm1 = 100
n_lstm2 = 100
n_hidden = 100
n_output = 4
model = tf.keras.models.Sequential([
tf.keras.layers.Input(shape=(train_window,4)), # (timestep, features)
tf.keras.layers.LSTM(n_lstm1, return_sequences = True),
tf.keras.layers.LSTM(n_lstm2),
tf.keras.layers.Dense(n_hidden),
tf.keras.layers.Dense(n_output),
])
model.summary()
model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=5e-4),
loss = 'mean_squared_error',
metrics = ['mse'])
In [ ]:
history = model.fit(train_data, train_label, epochs = 50, batch_size=32)
Evaluation¶
In [ ]:
# Prediction
preds = test_data[0].copy() # (window, 4)
for i in range(len(test_data)):
x = model.predict(np.expand_dims(preds[i:], axis = 0), verbose = 0) # (1, 4)
preds = np.concatenate((preds, x), axis = 0)
actual_predictions = scaler.inverse_transform(preds) # (1000,4)
In [ ]:
mae = mean_absolute_error(sol[550:], actual_predictions[50:])
mse = mean_squared_error(sol[550:], actual_predictions[50:])
print('Mean absolute error (MAE): {:.4f}'.format(mae))
print('Mean absolute error (MSE): {:.4f}'.format(mse))
In [ ]:
data = sol
data2 = np.concatenate((sol[:500], actual_predictions), axis = 0)
In [ ]:
def make_anim2(data1, data2):
def init():
line1.set_data([], [])
line2.set_data([], [])
line3.set_data([], [])
line4.set_data([], [])
line5.set_data([], [])
line6.set_data([], [])
line7.set_data([], [])
line8.set_data([], [])
line9.set_data([], [])
line10.set_data([], [])
time_string.set_text('')
return line3, line4, line5, line1, line2, line8, line9, line10, line6, line7, time_string
def animate(i):
trail1 = 6 # length of motion trail of weight 1
trail2 = 8 # length of motion trail of weight 2
dt = t[1] - t[0] # time step
# Update first set of lines
line1.set_data(xp1_1[i:max(1,i-trail1):-1], yp1_1[i:max(1,i-trail1):-1])
line2.set_data(xp2_1[i:max(1,i-trail2):-1], yp2_1[i:max(1,i-trail2):-1])
line3.set_data([xp1_1[i], xp2_1[i]], [yp1_1[i], yp2_1[i]])
line4.set_data([xp1_1[i], 0], [yp1_1[i], 0])
line5.set_data([0, 0], [0, 0])
# Update second set of lines
line6.set_data(xp1_2[i:max(1,i-trail1):-1], yp1_2[i:max(1,i-trail1):-1])
line7.set_data(xp2_2[i:max(1,i-trail2):-1], yp2_2[i:max(1,i-trail2):-1])
line8.set_data([xp1_2[i], xp2_2[i]], [yp1_2[i], yp2_2[i]])
line9.set_data([xp1_2[i], 0], [yp1_2[i], 0])
line10.set_data([0, 0], [0, 0])
time_string.set_text(time_template % (i * dt))
return line3, line4, line5, line1, line2, line8, line9, line10, line6, line7, time_string
fig = plt.figure(figsize = (10,8))
ax = plt.axes(xlim=(-L1-L2-0.5, L1+L2+0.5), ylim=(-2.5, 1.5))
# First set of lines (data1)
line1, = ax.plot([], [], 'o-', color='#d2eeff', markersize=12, markerfacecolor='#0077BE', lw=2,
markevery=10000, markeredgecolor='k')
line2, = ax.plot([], [], 'o-', color='#ffebd8', markersize=12, markerfacecolor='#f66338', lw=2,
markevery=10000, markeredgecolor='k', label = 'Ground Truth')
line3, = ax.plot([], [], color='k', linestyle='-', linewidth=2)
line4, = ax.plot([], [], color='k', linestyle='-', linewidth=2)
line5, = ax.plot([], [], 'o', color='k', markersize=10)
# Second set of lines (data2)
line6, = ax.plot([], [], 'o-', color='#AEEEEE', markersize=12, markerfacecolor='#5F9EA0', lw=2,
markevery=10000, markeredgecolor='k')
line7, = ax.plot([], [], 'o-', color='#FFC1C1', markersize=12, markerfacecolor='#DC143C', lw=2,
markevery=10000, markeredgecolor='k', label='Prediction')
line8, = ax.plot([], [], color='k', linestyle='-', linewidth=2)
line9, = ax.plot([], [], color='k', linestyle='-', linewidth=2)
line10, = ax.plot([], [], 'o', color='k', markersize=10)
time_template = 'Time = %.1f s'
time_string = ax.text(0.05, 0.9, '', transform=ax.transAxes, fontsize = 15)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
# Calculate coordinates for data1
u0_1 = data1[:, 0]
u1_1 = data1[:, 1]
u2_1 = data1[:, 2]
u3_1 = data1[:, 3]
xp1_1 = L1 * np.sin(u0_1)
yp1_1 = -L1 * np.cos(u0_1)
xp2_1 = xp1_1 + L2 * np.sin(u2_1)
yp2_1 = yp1_1 - L2 * np.cos(u2_1)
# Calculate coordinates for data2
u0_2 = data2[:, 0]
u1_2 = data2[:, 1]
u2_2 = data2[:, 2]
u3_2 = data2[:, 3]
xp1_2 = L1 * np.sin(u0_2)
yp1_2 = -L1 * np.cos(u0_2)
xp2_2 = xp1_2 + L2 * np.sin(u2_2)
yp2_2 = yp1_2 - L2 * np.cos(u2_2)
ax.legend(loc='upper right', fontsize = 15) # Add the legend to the plot
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=Nt, interval=1000 * (t[1] - t[0]) * 0.8, blit=True)
return anim
In [ ]:
anim_overlap = make_anim2(data,data2)
In [ ]:
HTML(anim_overlap.to_html5_video())
Out[ ]:
Lab 2: LSTM for Prediction¶
In [ ]:
theta1 = -60
theta2 = -30
state = np.radians([theta1, initial_velo_1, theta2, initial_velo_2])
sol = odeint(derivs, state, t)
In [ ]:
scaler = MinMaxScaler(feature_range=(-1, 1)) # -1 ~ 1 normalization
train_data_normalized = scaler.fit_transform(sol)
train_data_normalized = np.array(train_data_normalized)
data, label= create_inout_sequences(train_data_normalized, train_window)
train_data = np.array(data).reshape((-1,train_window,4))[:500]
test_data = np.array(data).reshape((-1,train_window,4))[500:]
train_label = np.array(label).reshape((-1,4))[:500]
test_label = np.array(label).reshape((-1,4))[500:]
In [ ]:
n_lstm1 = 100
n_lstm2 = 100
n_hidden = 100
n_output = 4
model = tf.keras.models.Sequential([
tf.keras.layers.Input(shape=(train_window,4)), # (timestep, features)
tf.keras.layers.LSTM(n_lstm1, return_sequences = True),
tf.keras.layers.LSTM(n_lstm2),
tf.keras.layers.Dense(n_hidden, activation = 'relu'),
tf.keras.layers.Dense(n_output),
])
model.summary()
model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=5e-4),
loss = 'mean_squared_error',
metrics = ['mse'])
In [ ]:
history = model.fit(train_data, train_label, epochs = 50, batch_size=32)
Evaluation¶
In [ ]:
# Prediction
preds = test_data[0].copy() # (window, 4)
for i in range(len(test_data)):
x = model.predict(np.expand_dims(preds[i:], axis = 0), verbose = 0) # (1, 4)
preds = np.concatenate((preds, x), axis = 0)
actual_predictions = scaler.inverse_transform(preds) # (1000,4)
In [ ]:
mae = mean_absolute_error(sol[550:], actual_predictions[50:])
mse = mean_squared_error(sol[550:], actual_predictions[50:])
print('Mean absolute error (MAE): {:.4f}'.format(mae))
print('Mean absolute error (MSE): {:.4f}'.format(mse))
In [ ]:
data = sol
data2 = np.concatenate((sol[:500], actual_predictions), axis = 0)
In [ ]:
anim_overlap = make_anim2(data,data2)
In [ ]:
HTML(anim_overlap.to_html5_video())
Out[ ]:
Lab 3: LSTM for Prediction¶
In [ ]:
theta1 = -30
theta2 = 30
state = np.radians([theta1, initial_velo_1, theta2, initial_velo_2])
sol = odeint(derivs, state, t)
In [ ]:
scaler = MinMaxScaler(feature_range=(-1, 1)) # -1 ~ 1 normalization
train_data_normalized = scaler.fit_transform(sol)
train_data_normalized = np.array(train_data_normalized)
data, label= create_inout_sequences(train_data_normalized, train_window)
train_data = np.array(data).reshape((-1,train_window,4))[:500]
test_data = np.array(data).reshape((-1,train_window,4))[500:]
train_label = np.array(label).reshape((-1,4))[:500]
test_label = np.array(label).reshape((-1,4))[500:]
In [ ]:
n_lstm1 = 100
n_lstm2 = 100
n_hidden = 100
n_output = 4
model = tf.keras.models.Sequential([
tf.keras.layers.Input(shape=(train_window,4)), # (timestep, features)
tf.keras.layers.LSTM(n_lstm1, return_sequences = True),
tf.keras.layers.LSTM(n_lstm2),
tf.keras.layers.Dense(n_hidden, activation = 'relu'),
tf.keras.layers.Dense(n_output),
])
model.summary()
model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=5e-4),
loss = 'mean_squared_error',
metrics = ['mse'])
In [ ]:
history = model.fit(train_data, train_label, epochs = 50, batch_size=32)
Evaluation¶
In [ ]:
# Prediction
preds = test_data[0].copy() # (window, 4)
for i in range(len(test_data)):
x = model.predict(np.expand_dims(preds[i:], axis = 0), verbose = 0) # (1, 4)
preds = np.concatenate((preds, x), axis = 0)
actual_predictions = scaler.inverse_transform(preds) # (1000,4)
In [ ]:
mae = mean_absolute_error(sol[550:], actual_predictions[50:])
mse = mean_squared_error(sol[550:], actual_predictions[50:])
print('Mean absolute error (MAE): {:.4f}'.format(mae))
print('Mean absolute error (MSE): {:.4f}'.format(mse))
In [ ]:
data = sol
data2 = np.concatenate((sol[:500], actual_predictions), axis = 0)
In [ ]:
anim_overlap = make_anim2(data,data2)
In [ ]:
HTML(anim_overlap.to_html5_video())
Out[ ]:
Vehicle Speed Estimation¶
Data download: - Vehicle video
Objective
- Develop a deep learning model to accurately estimate vehicle velocity.
- Output precise velocity predictions to enhance autonomous driving systems.
Limitations of Traditional Methods
- Sensor Dependency: Rely heavily on GPS, LIDAR, and RADAR, which can be expensive and prone to errors in certain conditions
- Speed Calculation using Wheel Rotation: Slippage, skidding, tire wear, etc.
In [ ]:
import random
from collections import OrderedDict
from pathlib import Path
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
from natsort import natsorted
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import *
from tensorflow.keras.applications import *
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.losses import MeanSquaredError
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.utils import Sequence
import cv2
import os
from pathlib import Path
In [ ]:
# Setting random seed
random.seed(hash("setting random seeds") % 2**32 - 1)
# Configuration dictionary
p = {
'project_name': 'vehicle-speed-estimation',
'batch_size': 32,
'w': 240,
'h': 320,
'model': 'efficientnetb0',
'split': 0.2,
'mean': .1,
'std': .5,
'divide_y': 10
}
# File paths for images and labels
IMAGES = './data/frames_test/'
LABELS = './data/train_vel_test.txt'
In [ ]:
def preprocess_image(image, w, h, mean, std):
# Resize
image = tf.image.resize(image, [w, h])
# Rescale
image = image / 255.0
# Normalize
image = (image - mean) / std
return image
class DS(tf.keras.utils.Sequence):
def __init__(self, images, labels, w, h, mean, std):
self.images = Path(images)
self.labels = [float(line.split()[0]) for line in open(labels).readlines()]
self.n_images = len(list(self.images.glob('*.png')))
self.indices = list(range(self.n_images - 1))
self.batch_size = p['batch_size']
self.w = w
self.h = h
self.mean = mean
self.std = std
def __len__(self):
return int(np.ceil(len(self.indices) / self.batch_size))
def __getitem__(self, idx):
batch_indices = self.indices[idx * self.batch_size:(idx + 1) * self.batch_size]
x = []
y = []
for index in batch_indices:
f1 = f'{self.images}/{str(index)}.png'
f2 = f'{self.images}/{str(index+1)}.png'
image1 = Image.open(f1)
image2 = Image.open(f2)
image1 = preprocess_image(np.array(image1), self.w, self.h, self.mean, self.std)
image2 = preprocess_image(np.array(image2), self.w, self.h, self.mean, self.std)
x.append(tf.concat((image1, image2), axis=-1))
y.append(self.labels[index] / p['divide_y'])
return tf.stack(x), tf.convert_to_tensor(y)
# Example usage of DS class
train_ds = DS(IMAGES, LABELS, p['w'], p['h'], p['mean'], p['std'])
valid_ds = DS(IMAGES, LABELS, p['w'], p['h'], p['mean'], p['std'])
test_ds = DS('./data/frames_test/', './data/train_vel_test.txt', p['w'], p['h'], p['mean'], p['std'])
In [ ]:
image_files = natsorted(os.listdir(IMAGES))
X = np.zeros((len(image_files)-1, p['w'], p['h'], 6))
y = np.array([float(line.split()[0]) for line in open(LABELS).readlines()])
for i in range(len(image_files)-1):
image1 = Image.open(IMAGES+image_files[i])
image2 = Image.open(IMAGES+image_files[i+1])
image1 = preprocess_image(np.array(image1), p['w'], p['h'], p['mean'], p['std'])
image2 = preprocess_image(np.array(image2), p['w'], p['h'], p['mean'], p['std'])
X[i] = tf.concat((image1, image2), axis = -1)
- Input: A sequence of 2 past frames of road scene images.
- Output: The predicted velocity of the vehicle.
In [ ]:
# Model input
input_layer = Input(shape=(p['w'], p['h'], 6))
x1 = Conv2D(filters = 32, kernel_size = (3,3), activation = 'relu', padding = 'SAME')(input_layer)
x1 = MaxPool2D((2,2))(x1)
x1 = Conv2D(filters = 64, kernel_size = (3,3), activation = 'relu', padding = 'SAME')(x1)
x1 = MaxPool2D((2,2))(x1)
x1 = Conv2D(filters = 128, kernel_size = (3,3), activation = 'relu', padding = 'SAME')(x1)
x1 = MaxPool2D((2,2))(x1)
x1 = Conv2D(filters = 256, kernel_size = (3,3), activation = 'relu', padding = 'SAME')(x1)
x1 = MaxPool2D((2,2))(x1)
x1 = Conv2D(filters = 512, kernel_size = (3,3), activation = 'relu', padding = 'SAME')(x1)
x1 = MaxPool2D((2,2))(x1)
x_f = Flatten()(x1)
x2 = Dense(128)(x_f)
encoded = Dense(6)(x2)
# Output layer
output_layer = Dense(1)(encoded)
# Build and compile the model
model = Model(inputs=input_layer, outputs=output_layer)
model.compile(optimizer=Adam(), loss=MeanSquaredError(), metrics=['mae'])
In [ ]:
model.summary()
In [ ]:
model.fit(train_ds, validation_data=valid_ds, epochs=30)
Out[ ]:
Evaluation¶
In [ ]:
image_files = natsorted(os.listdir('./data/frames_test/'))
X = np.zeros((len(image_files)-1, p['w'], p['h'], 6))
y = np.array([float(line.split()[0]) for line in open('./data/train_vel_test.txt').readlines()])
for i in range(len(image_files)-1):
image1 = Image.open('./data/frames_test/'+image_files[i])
image2 = Image.open('./data/frames_test/'+image_files[i+1])
image1 = preprocess_image(np.array(image1), p['w'], p['h'], p['mean'], p['std'])
image2 = preprocess_image(np.array(image2), p['w'], p['h'], p['mean'], p['std'])
X[i] = tf.concat((image1, image2), axis = -1)
pred = model.predict(X, verbose = 0)
plt.plot(pred*p['divide_y'], 'red')
plt.plot(y[:-1])
Out[ ]:
Visualization¶
In [ ]:
import os
from pathlib import Path
import cv2
import imageio
def annotate_frame_with_prediction(frame, prediction, ground_truth):
"""Annotate the frame with the prediction and ground truth values."""
font = cv2.FONT_HERSHEY_SIMPLEX
org_pred = (50, 50)
org_gt = (50, 100)
font_scale = 1
color_pred = (255, 0, 0)
color_gt = (0, 255, 0)
thickness = 2
annotated_frame = cv2.putText(frame, f'Prediction: {prediction:.2f} mph', org_pred, font, font_scale, color_pred, thickness, cv2.LINE_AA)
annotated_frame = cv2.putText(annotated_frame, f'Ground Truth: {ground_truth:.2f} mph', org_gt, font, font_scale, color_gt, thickness, cv2.LINE_AA)
return annotated_frame
def create_gif_with_predictions(input_dir, labels_file, output_gif, model):
"""Create a GIF from annotated frames with predictions and ground truth."""
frames = []
frame_files = sorted([f for f in os.listdir(input_dir) if f.endswith('.png')], key=lambda x: int(x.split('.')[0]))
labels = [float(line.split()[0]) for line in open(labels_file).readlines()]
for i in range(len(frame_files) - 1):
frame_file_1 = frame_files[i]
frame_file_2 = frame_files[i + 1]
frame_path_1 = os.path.join(input_dir, frame_file_1)
frame_path_2 = os.path.join(input_dir, frame_file_2)
frame_1 = cv2.imread(frame_path_1)
frame_2 = cv2.imread(frame_path_2)
# Preprocess the frames as needed for the model
frame_1_resized = preprocess_image(frame_1, p['w'], p['h'], p['mean'], p['std'])
frame_2_resized = preprocess_image(frame_2, p['w'], p['h'], p['mean'], p['std'])
# Combine frames along the channel dimension
combined_frames = tf.concat((frame_1_resized, frame_2_resized), axis=-1)
frame_input = tf.reshape(combined_frames, (1, p['w'], p['h'], 3*2))
# Make prediction
prediction = model.predict(frame_input)[0][0]
# Ground truth
ground_truth = labels[i] / p['divide_y']
# Annotate the first frame with prediction and ground truth
annotated_frame = annotate_frame_with_prediction(frame_1, prediction*10, ground_truth*10)
resized_annotated_frame = cv2.resize(annotated_frame, (p['w'], p['h']))
# Convert BGR (OpenCV format) to RGB (imageio format)
annotated_frame_rgb = cv2.cvtColor(resized_annotated_frame, cv2.COLOR_BGR2RGB)
frames.append(annotated_frame_rgb)
# Save frames as a GIF
imageio.mimsave(output_gif, frames, fps=20)
# Example usage:
input_dir = "./data/frames_test"
labels_file = "./data/train_vel_test.txt"
output_gif = "./data/output.gif"
create_gif_with_predictions(input_dir, labels_file, output_gif, model)