1777 lines
78 KiB
Plaintext
1777 lines
78 KiB
Plaintext
[1]. Tool: Gemini
|
|
Prompt: Could you please set up my RL training environment to load the URDF file I created for my maglev pod and
|
|
put a track above it? This should all be done in PyBullet simulation, using the maglev_predictor and
|
|
mag_lev_coil files I have provided. The agent should have access to an observation space with 4 gap
|
|
heights and 4 gap height velocities from the sim, and should be able to output 4 PWM values as actions.
|
|
Output:
|
|
# --- 1. Define Action & Observation Spaces ---
|
|
# Action: 4 PWM duty cycles between -1 and 1 (4 independent coils)
|
|
# [front_left, front_right, back_left, back_right] corresponding to +Y and -Y ends of each U-yoke
|
|
self.action_space = spaces.Box(low=-1, high=1, shape=(4,), dtype=np.float32)
|
|
|
|
# Observation: 4 normalized noisy sensor gap heights + 4 normalized velocities
|
|
# Gaps normalized by 0.030m, velocities by 0.1 m/s
|
|
self.observation_space = spaces.Box(low=-5.0, high=5.0, shape=(8,), dtype=np.float32)
|
|
|
|
# --- 2. Setup Physics & Actuators ---
|
|
self.dt = 1./240. # PyBullet default timestep
|
|
self.coil_front_L = MagLevCoil(1.1, 0.0025, 12, 10.2) # From your specs
|
|
self.coil_front_R = MagLevCoil(1.1, 0.0025, 12, 10.2)
|
|
self.coil_back_L = MagLevCoil(1.1, 0.0025, 12, 10.2)
|
|
self.coil_back_R = MagLevCoil(1.1, 0.0025, 12, 10.2)
|
|
|
|
# Sensor noise parameters
|
|
self.sensor_noise_std = 0.0001 # 0.1mm standard deviation
|
|
|
|
# Normalization constants for observations
|
|
self.gap_scale = 0.015 # Normalize gaps by +-15mm max expected deviation from middle
|
|
self.velocity_scale = 0.1 # Normalize velocities by 0.1 m/s max expected velocity
|
|
|
|
# Maglev force/torque predictor
|
|
self.predictor = MaglevPredictor()
|
|
|
|
# Connect to PyBullet (DIRECT is faster for training, GUI for debugging)
|
|
self.client = p.connect(p.GUI if use_gui else p.DIRECT)
|
|
p.setAdditionalSearchPath(pybullet_data.getDataPath())
|
|
|
|
# Store references
|
|
self.trackId = None
|
|
self.podId = None
|
|
self.collision_local_positions = []
|
|
self.yoke_indices = [] # For force application
|
|
self.yoke_labels = []
|
|
self.sensor_indices = [] # For gap height measurement
|
|
self.sensor_labels = []
|
|
|
|
# For velocity calculation
|
|
self.prev_sensor_gaps = None
|
|
|
|
|
|
def reset(self, seed=None, options=None):
|
|
# Reset PyBullet simulation
|
|
p.resetSimulation(physicsClientId=self.client)
|
|
p.setGravity(0, 0, -9.81, physicsClientId=self.client)
|
|
p.setTimeStep(self.dt, physicsClientId=self.client)
|
|
|
|
# Create the maglev track (inverted system - track above, pod hangs below)
|
|
# Track bottom surface at Z=0
|
|
track_collision = p.createCollisionShape(
|
|
shapeType=p.GEOM_BOX,
|
|
halfExtents=[1.0, 0.2, 0.010],
|
|
physicsClientId=self.client
|
|
)
|
|
track_visual = p.createVisualShape(
|
|
shapeType=p.GEOM_BOX,
|
|
halfExtents=[1.0, 0.2, 0.010],
|
|
rgbaColor=[0.3, 0.3, 0.3, 0.8],
|
|
physicsClientId=self.client
|
|
)
|
|
self.trackId = p.createMultiBody(
|
|
baseMass=0, # Static
|
|
baseCollisionShapeIndex=track_collision,
|
|
baseVisualShapeIndex=track_visual,
|
|
basePosition=[0, 0, 0.010], # Track center at Z=10mm, bottom at Z=0
|
|
physicsClientId=self.client
|
|
)
|
|
p.changeDynamics(self.trackId, -1,
|
|
lateralFriction=0.3,
|
|
restitution=0.1,
|
|
physicsClientId=self.client)
|
|
|
|
start_pos = [0, 0, start_z]
|
|
start_orientation = p.getQuaternionFromEuler([0, 0, 0])
|
|
self.podId = p.loadURDF(urdf_path, start_pos, start_orientation, physicsClientId=self.client)
|
|
|
|
[2]. Tool: Gemini
|
|
Prompt:
|
|
Output:
|
|
# Parse collision shapes to identify yokes and sensors
|
|
collision_shapes = p.getCollisionShapeData(self.podId, -1, physicsClientId=self.client)
|
|
self.collision_local_positions = []
|
|
self.yoke_indices = []
|
|
self.yoke_labels = []
|
|
self.sensor_indices = []
|
|
self.sensor_labels = []
|
|
|
|
# Expected heights for detection
|
|
expected_yoke_sensor_z = 0.08585 # Yokes and sensors always at this height
|
|
expected_bolt_z = 0.08585 + self.initial_gap_mm / 1000.0 # Bolts at gap-dependent height
|
|
|
|
for i, shape in enumerate(collision_shapes):
|
|
shape_type = shape[2]
|
|
local_pos = shape[5]
|
|
self.collision_local_positions.append(local_pos)
|
|
|
|
# Check if at sensor/yoke height (Z ≈ 0.08585m) - NOT bolts
|
|
if abs(local_pos[2] - expected_yoke_sensor_z) < 0.001:
|
|
if shape_type == p.GEOM_BOX:
|
|
# Yokes are BOX type at the four corners (size 0.0254)
|
|
self.yoke_indices.append(i)
|
|
x_pos = "Front" if local_pos[0] > 0 else "Back"
|
|
y_pos = "Left" if local_pos[1] > 0 else "Right"
|
|
self.yoke_labels.append(f"{x_pos}_{y_pos}")
|
|
elif shape_type == p.GEOM_CYLINDER or shape_type == p.GEOM_MESH:
|
|
# Sensors: distinguish by position pattern
|
|
if abs(local_pos[0]) < 0.06 or abs(local_pos[1]) < 0.02:
|
|
self.sensor_indices.append(i)
|
|
if abs(local_pos[0]) < 0.001: # Center sensors (X ≈ 0)
|
|
label = "Center_Right" if local_pos[1] > 0 else "Center_Left"
|
|
else: # Front/back sensors (Y ≈ 0)
|
|
label = "Front" if local_pos[0] > 0 else "Back"
|
|
self.sensor_labels.append(label)
|
|
|
|
[3]. Tool: Claude Code
|
|
Prompt: Given the maglev_predictor file, generate code to find the equilibrium height for a pod of a given mass.
|
|
This should simply be the gap height at which 0 current in either coil and zero roll (all times 2) creates a
|
|
force equivalent to the gravitational force on the 5.8kg pod.
|
|
|
|
Output:
|
|
"""
|
|
Find equilibrium gap height for magnetic levitation system.
|
|
|
|
Given:
|
|
- Pod mass: 5.8 kg
|
|
- Required force: 5.8 * 9.81 = 56.898 N
|
|
- All currents: 0 A
|
|
- Roll angle: 0 degrees
|
|
|
|
Find: Gap height (mm) that produces this force
|
|
"""
|
|
|
|
import numpy as np
|
|
from maglev_predictor import MaglevPredictor
|
|
|
|
# Initialize predictor
|
|
predictor = MaglevPredictor()
|
|
|
|
# Target force
|
|
target_force = 5.8 * 9.81 # 56.898 N
|
|
|
|
print("=" * 70)
|
|
print("EQUILIBRIUM GAP HEIGHT FINDER")
|
|
print("=" * 70)
|
|
print(f"Target Force: {target_force:.3f} N (for 5.8 kg pod)")
|
|
print(f"Parameters: currL = 0 A, currR = 0 A, roll = 0°")
|
|
print()
|
|
|
|
# Define objective function
|
|
def force_error(gap_height):
|
|
"""Calculate difference between predicted force and target force."""
|
|
force, _ = predictor.predict(currL=0, currR=0, roll=0, gap_height=gap_height)
|
|
return force - target_force
|
|
|
|
# First, let's scan across gap heights to understand the behavior
|
|
print("Scanning gap heights from 5 to 30 mm...")
|
|
print("-" * 70)
|
|
gap_range = np.linspace(5, 30, 26)
|
|
forces = []
|
|
|
|
for gap in gap_range:
|
|
force, torque = predictor.predict(currL=0, currR=0, roll=0, gap_height=gap)
|
|
forces.append(force)
|
|
if abs(force - target_force) < 5: # Close to target
|
|
print(f"Gap: {gap:5.1f} mm → Force: {force:7.3f} N ← CLOSE!")
|
|
else:
|
|
print(f"Gap: {gap:5.1f} mm → Force: {force:7.3f} N")
|
|
|
|
forces = np.array(forces)
|
|
|
|
print()
|
|
print("-" * 70)
|
|
|
|
# Check if target force is within the range
|
|
min_force = forces.min()
|
|
max_force = forces.max()
|
|
|
|
if target_force < min_force or target_force > max_force:
|
|
print(f"⚠ WARNING: Target force {target_force:.3f} N is outside the range!")
|
|
print(f" Force range: {min_force:.3f} N to {max_force:.3f} N")
|
|
print(f" Cannot achieve equilibrium with 0A currents.")
|
|
else:
|
|
# Find the gap height using root finding
|
|
# Use brentq for robust bracketing (requires sign change)
|
|
|
|
# Find bracketing interval
|
|
idx_above = np.where(forces > target_force)[0]
|
|
idx_below = np.where(forces < target_force)[0]
|
|
|
|
if len(idx_above) > 0 and len(idx_below) > 0:
|
|
# Find the transition point
|
|
gap_low = gap_range[idx_below[-1]]
|
|
gap_high = gap_range[idx_above[0]]
|
|
|
|
print(f"✓ Target force is achievable!")
|
|
print(f" Bracketing interval: {gap_low:.1f} mm to {gap_high:.1f} mm")
|
|
print()
|
|
|
|
# Use simple bisection method for accurate root finding
|
|
tol = 1e-6
|
|
while (gap_high - gap_low) > tol:
|
|
gap_mid = (gap_low + gap_high) / 2
|
|
force_mid, _ = predictor.predict(currL=0, currR=0, roll=0, gap_height=gap_mid)
|
|
|
|
if force_mid > target_force:
|
|
gap_low = gap_mid
|
|
else:
|
|
gap_high = gap_mid
|
|
|
|
equilibrium_gap = (gap_low + gap_high) / 2
|
|
|
|
# Verify the result
|
|
final_force, final_torque = predictor.predict(
|
|
currL=0, currR=0, roll=0, gap_height=equilibrium_gap
|
|
)
|
|
|
|
print("=" * 70)
|
|
print("EQUILIBRIUM FOUND!")
|
|
print("=" * 70)
|
|
print(f"Equilibrium Gap Height: {equilibrium_gap:.6f} mm")
|
|
print(f"Predicted Force: {final_force:.6f} N")
|
|
print(f"Target Force: {target_force:.6f} N")
|
|
print(f"Error: {abs(final_force - target_force):.9f} N")
|
|
print(f"Torque at equilibrium: {final_torque:.6f} mN·m")
|
|
print()
|
|
print(f"✓ Pod will levitate at {equilibrium_gap:.3f} mm gap height")
|
|
print(f" with no current applied (permanent magnets only)")
|
|
print("=" * 70)
|
|
else:
|
|
print("⚠ Could not find bracketing interval for bisection.")
|
|
print(" Target force may not be achievable in the scanned range.")
|
|
|
|
[4]. Tool: Claude Code
|
|
Prompt:
|
|
How can I fit a polynomial to a system with 4 inputs and two outputs? I need to map the values
|
|
in the provided CSV file (currL, currR, roll, and inverse Gap) to Force and Torque readings.
|
|
|
|
Output:
|
|
# 1. Load your Ansys CSV - using invGap for modeling
|
|
X = magDf[['currL [A]', 'currR [A]', 'rollDeg [deg]', 'invGap']]
|
|
y = magDf[['YokeForce.Force_z [newton]', 'YokeTorque.Torque [mNewtonMeter]']]
|
|
|
|
# 2. Create Features (e.g. z^2, z^3, z*IL, etc.)
|
|
poly = PolynomialFeatures(degree=3)
|
|
X_poly = poly.fit_transform(X)
|
|
|
|
# 3. Fit
|
|
model = LinearRegression()
|
|
model.fit(X_poly, y)
|
|
|
|
# 4. Extract Equation
|
|
print("Force Coeffs:", model.coef_[0])
|
|
print("Torque Coeffs:", model.coef_[1])
|
|
|
|
[5]. Tool: Claude Code
|
|
Prompt:
|
|
Can you generate code for comparing a polynomial fit on the entire dataset (the actual
|
|
polynomial curve, not just the points that are provided by Ansys) with the ansys points
|
|
themselves for a single segment with the currents and roll degree held constant and the
|
|
Gap Height on X-axis and force on y-axis? Use the given code for reference.
|
|
|
|
Output:
|
|
# Train best degree polynomial on full dataset
|
|
poly_best = PolynomialFeatures(degree=best_degree)
|
|
X_poly_best = poly_best.fit_transform(X)
|
|
model_best = LinearRegression()
|
|
model_best.fit(X_poly_best, y)
|
|
|
|
print(f"Using polynomial degree: {best_degree}")
|
|
|
|
# Select a specific operating condition to visualize
|
|
# Using the same condition as earlier: -15A, -15A, 0° roll
|
|
test_condition = magDf.iloc[0:13].copy()
|
|
currL_val = test_condition['currL [A]'].iloc[0]
|
|
currR_val = test_condition['currR [A]'].iloc[0]
|
|
roll_val = test_condition['rollDeg [deg]'].iloc[0]
|
|
|
|
# Create very fine grid for gap heights to see polynomial behavior between points
|
|
gap_very_fine = np.linspace(5, 30, 500) # 500 points for smooth curve
|
|
X_fine = pd.DataFrame({
|
|
'currL [A]': [currL_val] * 500,
|
|
'currR [A]': [currR_val] * 500,
|
|
'rollDeg [deg]': [roll_val] * 500,
|
|
'invGap': 1/gap_very_fine # Use inverse gap height for modeling
|
|
})
|
|
X_fine_poly = poly_best.transform(X_fine)
|
|
y_fine_pred = model_best.predict(X_fine_poly)
|
|
|
|
# Extract actual Ansys data for this condition
|
|
gap_actual = test_condition['GapHeight [mm]'].values
|
|
force_actual = test_condition['YokeForce.Force_z [newton]'].values
|
|
|
|
# Plot
|
|
fig, ax = plt.subplots(1, 1, figsize=(12, 7))
|
|
|
|
ax.plot(gap_very_fine, y_fine_pred[:, 0], '-', linewidth=2.5,
|
|
label=f'Degree {best_degree} Polynomial Function', color='#2ca02c', alpha=0.8)
|
|
|
|
# Plot actual Ansys data points
|
|
ax.scatter(gap_actual, force_actual, s=120, marker='o',
|
|
color='#d62728', edgecolors='black', linewidths=1.5,
|
|
label='Ansys Simulation Data', zorder=5)
|
|
|
|
ax.set_ylabel('Force (N)', fontsize=13)
|
|
ax.set_title(f'Degree {best_degree} Polynomial vs Ansys Data (currL={currL_val}A, currR={currR_val}A, roll={roll_val}°)',
|
|
fontsize=14, fontweight='bold')
|
|
ax.legend(fontsize=12, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
|
|
# Add vertical lines at data points to highlight gaps
|
|
for gap_point in gap_actual:
|
|
ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
# Check for excessive oscillations by looking at derivative
|
|
print("Checking for overfitting indicators:")
|
|
print(f"Number of actual data points in this condition: {len(gap_actual)}")
|
|
|
|
print(f"Gap spacing: {np.diff(gap_actual).mean():.2f} mm average")
|
|
print(f"Polynomial curve resolution: 500 points over {gap_very_fine.max() - gap_very_fine.min():.1f} mm range")
|
|
|
|
[6]. Tool: Claude Code
|
|
Prompt:
|
|
Now can we do four different trials with similar graphs but different combinations of the starting conditions?
|
|
Output:
|
|
# Select 4 different conditions to visualize
|
|
# Each condition has 13 gap height measurements (one complete parameter set)
|
|
condition_indices = [
|
|
(0, 13, '-15A, -15A, 0°'), # First condition
|
|
(13, 26, '-15A, -15A, 0.5°'), # Same currents, different roll
|
|
(91, 104, '-15A, -10A, 0°'), # Different right current
|
|
(182, 195, '-15A, -10A, 1°') # Different right current and roll
|
|
]
|
|
|
|
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
|
|
axes = axes.flatten()
|
|
|
|
for idx, (start, end, label) in enumerate(condition_indices):
|
|
ax = axes[idx]
|
|
|
|
# Get condition data
|
|
condition_data = magDf.iloc[start:end]
|
|
currL = condition_data['currL [A]'].iloc[0]
|
|
currR = condition_data['currR [A]'].iloc[0]
|
|
roll = condition_data['rollDeg [deg]'].iloc[0]
|
|
|
|
# Create fine grid
|
|
gap_fine = np.linspace(5, 30, 500)
|
|
X_condition_fine = pd.DataFrame({
|
|
'currL [A]': [currL] * 500,
|
|
'currR [A]': [currR] * 500,
|
|
'rollDeg [deg]': [roll] * 500,
|
|
'invGap': 1/gap_fine # Use inverse gap height for modeling
|
|
})
|
|
|
|
# Predict with best degree polynomial
|
|
X_condition_poly = poly_best.transform(X_condition_fine)
|
|
y_condition_pred = model_best.predict(X_condition_poly)
|
|
|
|
# Plot polynomial curve
|
|
ax.plot(gap_fine, y_condition_pred[:, 0], '-', linewidth=2.5,
|
|
color='#2ca02c', alpha=0.8, label=f'Degree {best_degree} Polynomial')
|
|
|
|
# Plot actual data
|
|
ax.scatter(condition_data['GapHeight [mm]'],
|
|
condition_data['YokeForce.Force_z [newton]'],
|
|
s=100, marker='o', color='#d62728',
|
|
edgecolors='black', linewidths=1.5,
|
|
label='Ansys Data', zorder=5)
|
|
|
|
# Formatting
|
|
ax.set_xlabel('Gap Height (mm)', fontsize=11)
|
|
ax.set_ylabel('Force (N)', fontsize=11)
|
|
ax.set_title(f'currL={currL:.0f}A, currR={currR:.0f}A, roll={roll:.0f}°',
|
|
fontsize=12, fontweight='bold')
|
|
ax.legend(fontsize=10, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
|
|
# Add vertical lines at data points
|
|
for gap_point in condition_data['GapHeight [mm]'].values:
|
|
ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)
|
|
|
|
plt.suptitle(f'Degree {best_degree} Polynomial: Multiple Operating Conditions',
|
|
fontsize=15, fontweight='bold', y=1.00)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print("Visual inspection guide:")
|
|
print("✓ Look for smooth curves between data points (good)")
|
|
print("✗ Look for wild oscillations or wiggles between points (overfitting)")
|
|
print("✓ Polynomial should interpolate naturally without excessive curvature")
|
|
|
|
[7]. Tool: Claude Code
|
|
Prompt:
|
|
Let's now find the worst performing segment of 13 gap heights to show the worst possible
|
|
performance of our polynomial fit in terms of the R^2 score
|
|
|
|
Output:
|
|
# Calculate R² score for each 13-row segment (each parameter set)
|
|
num_segments = len(magDf) // 13
|
|
segment_scores = []
|
|
segment_info = []
|
|
|
|
for i in range(num_segments):
|
|
start_idx = i * 13
|
|
end_idx = start_idx + 13
|
|
|
|
# Get actual and predicted values for this segment
|
|
segment_actual = y.iloc[start_idx:end_idx, 0].values # Force only
|
|
segment_pred = model_best.predict(poly_best.transform(X.iloc[start_idx:end_idx]))[:, 0]
|
|
|
|
# Calculate R² for this segment
|
|
segment_r2 = r2_score(segment_actual, segment_pred)
|
|
segment_scores.append(segment_r2)
|
|
|
|
# Store segment info
|
|
segment_data = magDf.iloc[start_idx:end_idx]
|
|
currL = segment_data['currL [A]'].iloc[0]
|
|
currR = segment_data['currR [A]'].iloc[0]
|
|
roll = segment_data['rollDeg [deg]'].iloc[0]
|
|
segment_info.append({
|
|
'start': start_idx,
|
|
'end': end_idx,
|
|
'currL': currL,
|
|
'currR': currR,
|
|
'roll': roll,
|
|
'r2': segment_r2
|
|
})
|
|
|
|
# Find worst performing segment
|
|
worst_idx = np.argmin(segment_scores)
|
|
worst_segment = segment_info[worst_idx]
|
|
|
|
print(f"Total segments analyzed: {num_segments}")
|
|
print(f"\nWorst performing segment:")
|
|
print(f" Rows: {worst_segment['start']} to {worst_segment['end']}")
|
|
print(f" Parameters: currL={worst_segment['currL']:.0f}A, currR={worst_segment['currR']:.0f}A, roll={worst_segment['roll']:.1f}°")
|
|
print(f" R² Score: {worst_segment['r2']:.6f}")
|
|
print(f"\nBest performing segment R²: {max(segment_scores):.6f}")
|
|
print(f"Mean R² across all segments: {np.mean(segment_scores):.6f}")
|
|
print(f"Median R² across all segments: {np.median(segment_scores):.6f}")
|
|
|
|
[8]. Tool: Claude Code
|
|
Prompt:
|
|
Let's also visualize the worst segment so the evidently good performance is visible to a viewer.
|
|
|
|
Output:
|
|
# Get worst segment data
|
|
worst_data = magDf.iloc[worst_segment['start']:worst_segment['end']]
|
|
worst_currL = worst_segment['currL']
|
|
worst_currR = worst_segment['currR']
|
|
worst_roll = worst_segment['roll']
|
|
worst_r2 = worst_segment['r2']
|
|
|
|
# Create fine grid for this worst segment
|
|
gap_fine_worst = np.linspace(5, 30, 500)
|
|
X_worst_fine = pd.DataFrame({
|
|
'currL [A]': [worst_currL] * 500,
|
|
'currR [A]': [worst_currR] * 500,
|
|
'rollDeg [deg]': [worst_roll] * 500,
|
|
'invGap': 1/gap_fine_worst
|
|
})
|
|
|
|
# Get predictions
|
|
X_worst_poly = poly_best.transform(X_worst_fine)
|
|
y_worst_pred = model_best.predict(X_worst_poly)
|
|
|
|
# Extract actual data
|
|
gap_worst_actual = worst_data['GapHeight [mm]'].values
|
|
force_worst_actual = worst_data['YokeForce.Force_z [newton]'].values
|
|
|
|
# Create the plot
|
|
fig, ax = plt.subplots(1, 1, figsize=(12, 7))
|
|
|
|
# Plot polynomial curve
|
|
ax.plot(gap_fine_worst, y_worst_pred[:, 0], '-', linewidth=2.5,
|
|
label=f'Degree {best_degree} Polynomial', color='#2ca02c', alpha=0.8)
|
|
|
|
# Plot actual data points
|
|
ax.scatter(gap_worst_actual, force_worst_actual, s=120, marker='o',
|
|
color='#d62728', edgecolors='black', linewidths=1.5,
|
|
label='Ansys Data', zorder=5)
|
|
|
|
# Formatting
|
|
ax.set_xlabel('Gap Height (mm)', fontsize=13)
|
|
ax.set_ylabel('Force (N)', fontsize=13)
|
|
ax.set_title(f'Worst Case Performance (R² = {worst_r2:.6f}) - currL={worst_currL:.0f}A, currR={worst_currR:.0f}A, roll={worst_roll:.1f}°',
|
|
fontsize=14, fontweight='bold')
|
|
ax.legend(fontsize=12, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
|
|
# Add vertical lines at data points
|
|
for gap_point in gap_worst_actual:
|
|
ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print(f"Even in the worst case, the model achieves R² = {worst_r2:.6f}")
|
|
print(f"This demonstrates excellent interpolation performance across all operating conditions!")
|
|
|
|
[9]. Tool: Claude Code
|
|
Prompt:
|
|
Now let's do something similar for the torque, showing how our polynomial fits to the torque vs
|
|
roll trend in ansys data while keeping othe variables constant. Show 4 trials like before, and
|
|
generate a continuous plot, like before.
|
|
|
|
Output:
|
|
# Select 4 different gap heights to visualize
|
|
# Using currL = -15A, currR = -15A configuration
|
|
gap_heights_to_plot = [8, 10, 15, 21] # mm
|
|
|
|
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
|
|
axes = axes.flatten()
|
|
|
|
for idx, gap_height in enumerate(gap_heights_to_plot):
|
|
ax = axes[idx]
|
|
|
|
# Fixed current configuration
|
|
currL = -15
|
|
currR = -15
|
|
|
|
# Get actual data for this gap height
|
|
gap_data = magDf[(magDf['GapHeight [mm]'] == gap_height) &
|
|
(magDf['currL [A]'] == currL) &
|
|
(magDf['currR [A]'] == currR)]
|
|
|
|
# Create fine grid for smooth curve across roll angles
|
|
roll_fine = np.linspace(-5, 5.0, 500)
|
|
X_condition_fine = pd.DataFrame({
|
|
'currL [A]': [currL] * 500,
|
|
'currR [A]': [currR] * 500,
|
|
'rollDeg [deg]': roll_fine,
|
|
'invGap': [1/gap_height] * 500
|
|
})
|
|
|
|
# Predict with best degree polynomial
|
|
X_condition_poly = poly_best.transform(X_condition_fine)
|
|
y_condition_pred = model_best.predict(X_condition_poly)
|
|
|
|
# Plot polynomial curve for TORQUE (index 1)
|
|
ax.plot(roll_fine, y_condition_pred[:, 1], '-', linewidth=2.5,
|
|
color='#ff7f0e', alpha=0.8, label=f'Degree {best_degree} Polynomial')
|
|
|
|
# Plot actual torque data
|
|
ax.scatter(gap_data['rollDeg [deg]'],
|
|
gap_data['YokeTorque.Torque [mNewtonMeter]'],
|
|
s=100, marker='o', color='#9467bd',
|
|
edgecolors='black', linewidths=1.5,
|
|
label='Ansys Data', zorder=5)
|
|
|
|
# Formatting
|
|
ax.set_xlabel('Roll Angle (deg)', fontsize=11)
|
|
ax.set_ylabel('Torque (mN·m)', fontsize=11)
|
|
ax.set_title(f'Gap Height = {gap_height} mm (currL={currL:.0f}A, currR={currR:.0f}A)',
|
|
fontsize=12, fontweight='bold')
|
|
ax.legend(fontsize=10, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
|
|
# Add vertical lines at actual roll angle data points
|
|
for roll_point in gap_data['rollDeg [deg]'].unique():
|
|
ax.axvline(roll_point, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)
|
|
|
|
plt.suptitle(f'Torque vs Roll Angle at Different Gap Heights (Degree {best_degree} Polynomial)',
|
|
fontsize=15, fontweight='bold', y=1.00)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print(f"Showing torque variation with roll angle at {len(gap_heights_to_plot)} different gap heights.")
|
|
print("Gap height affects the magnitude of torque but the relationship with roll angle remains consistent.")
|
|
|
|
[10]. Tool: Claude Code
|
|
Prompt:
|
|
Since torque also varies with current imbalance between the coils, let's visualize the
|
|
response for that as well. Follow the same visualization pattern as before.
|
|
|
|
Output:
|
|
# Select 4 different gap heights to visualize
|
|
# Using roll = 0.5 degrees (small but non-zero torque)
|
|
gap_heights_to_plot = [8, 10, 15, 21] # mm
|
|
roll_target = 0.5
|
|
|
|
# Use a fixed reference current for currL
|
|
currL_ref = -15
|
|
|
|
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
|
|
axes = axes.flatten()
|
|
|
|
for idx, gap_height in enumerate(gap_heights_to_plot):
|
|
ax = axes[idx]
|
|
|
|
# Get actual data for this gap height, roll angle, AND fixed currL
|
|
gap_roll_data = magDf[(magDf['GapHeight [mm]'] == gap_height) &
|
|
(magDf['rollDeg [deg]'] == roll_target) &
|
|
(magDf['currL [A]'] == currL_ref)]
|
|
|
|
# Calculate imbalance for actual data points (all have same currL now)
|
|
actual_imbalances = []
|
|
actual_torques = []
|
|
for _, row in gap_roll_data.iterrows():
|
|
imbalance = row['currR [A]'] - row['currL [A]']
|
|
torque = row['YokeTorque.Torque [mNewtonMeter]']
|
|
actual_imbalances.append(imbalance)
|
|
actual_torques.append(torque)
|
|
|
|
# Create fine grid for smooth curve across current imbalances
|
|
# Range from -10A to +10A imbalance (currR - currL)
|
|
imbalance_fine = np.linspace(0, 35, 500)
|
|
|
|
# Vary currR while keeping currL fixed
|
|
currR_fine = currL_ref + imbalance_fine
|
|
|
|
X_condition_fine = pd.DataFrame({
|
|
'currL [A]': [currL_ref] * 500,
|
|
'currR [A]': currR_fine,
|
|
'rollDeg [deg]': [roll_target] * 500,
|
|
'invGap': [1/gap_height] * 500
|
|
})
|
|
|
|
# Predict with best degree polynomial
|
|
X_condition_poly = poly_best.transform(X_condition_fine)
|
|
y_condition_pred = model_best.predict(X_condition_poly)
|
|
|
|
# Plot polynomial curve for TORQUE (index 1)
|
|
ax.plot(imbalance_fine, y_condition_pred[:, 1], '-', linewidth=2.5,
|
|
color='#ff7f0e', alpha=0.8, label=f'Degree {best_degree} Polynomial')
|
|
|
|
# Plot actual torque data (now only one point per imbalance value)
|
|
ax.scatter(actual_imbalances, actual_torques,
|
|
s=100, marker='o', color='#9467bd',
|
|
edgecolors='black', linewidths=1.5,
|
|
label='Ansys Data', zorder=5)
|
|
|
|
# Formatting
|
|
ax.set_xlabel('Current Imbalance (currR - currL) [A]', fontsize=11)
|
|
ax.set_ylabel('Torque (mN·m)', fontsize=11)
|
|
ax.set_title(f'Gap Height = {gap_height} mm (currL = {currL_ref}A, roll = {roll_target}°)',
|
|
fontsize=12, fontweight='bold')
|
|
ax.legend(fontsize=10, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
ax.axhline(0, color='black', linestyle='-', linewidth=0.5, alpha=0.3)
|
|
ax.axvline(0, color='black', linestyle='-', linewidth=0.5, alpha=0.3)
|
|
|
|
# Add vertical lines at actual imbalance data points
|
|
for imb in set(actual_imbalances):
|
|
ax.axvline(imb, color='gray', linestyle=':', alpha=0.2, linewidth=0.8)
|
|
|
|
plt.suptitle(f'Torque vs Current Imbalance at Different Gap Heights (currL = {currL_ref}A, Roll = {roll_target}°)',
|
|
fontsize=15, fontweight='bold', y=1.00)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print(f"Showing torque variation with current imbalance at {len(gap_heights_to_plot)} different gap heights.")
|
|
print(f"All data uses currL = {currL_ref}A with varying currR to create different imbalances.")
|
|
print(f"Current imbalance (currR - currL) affects torque magnitude and direction.")
|
|
print(f"Positive imbalance (currR > currL) and negative imbalance (currR < currL) produce opposite torques.")
|
|
|
|
[11]. Tool: Claude Code
|
|
Prompt:
|
|
Let's now visualize the worst torque prediction for a given 13-gap-height segment.
|
|
What is our worst possible performance? Show in a similar way as with force.
|
|
|
|
Output:
|
|
# Calculate R² score for torque in each 13-row segment
|
|
torque_segment_scores = []
|
|
torque_segment_info = []
|
|
|
|
for i in range(num_segments):
|
|
start_idx = i * 13
|
|
end_idx = start_idx + 13
|
|
|
|
# Get actual and predicted torque values for this segment
|
|
segment_actual_torque = y.iloc[start_idx:end_idx, 1].values # Torque column
|
|
segment_pred_torque = model_best.predict(poly_best.transform(X.iloc[start_idx:end_idx]))[:, 1]
|
|
|
|
# Calculate R² for this segment
|
|
segment_r2_torque = r2_score(segment_actual_torque, segment_pred_torque)
|
|
torque_segment_scores.append(segment_r2_torque)
|
|
|
|
# Store segment info
|
|
segment_data = magDf.iloc[start_idx:end_idx]
|
|
currL = segment_data['currL [A]'].iloc[0]
|
|
currR = segment_data['currR [A]'].iloc[0]
|
|
roll = segment_data['rollDeg [deg]'].iloc[0]
|
|
torque_segment_info.append({
|
|
'start': start_idx,
|
|
'end': end_idx,
|
|
'currL': currL,
|
|
'currR': currR,
|
|
'roll': roll,
|
|
'r2': segment_r2_torque
|
|
})
|
|
|
|
# Find worst and best performing segments for torque
|
|
worst_torque_idx = np.argmin(torque_segment_scores)
|
|
best_torque_idx = np.argmax(torque_segment_scores)
|
|
worst_torque_segment = torque_segment_info[worst_torque_idx]
|
|
best_torque_segment = torque_segment_info[best_torque_idx]
|
|
|
|
print(f"Torque Prediction Analysis (across {num_segments} segments):")
|
|
print(f"\nBest R² Score: {max(torque_segment_scores):.6f}")
|
|
print(f" Parameters: currL={best_torque_segment['currL']:.0f}A, currR={best_torque_segment['currR']:.0f}A, roll={best_torque_segment['roll']:.1f}°")
|
|
print(f"\nWorst R² Score: {min(torque_segment_scores):.6f}")
|
|
print(f" Parameters: currL={worst_torque_segment['currL']:.0f}A, currR={worst_torque_segment['currR']:.0f}A, roll={worst_torque_segment['roll']:.1f}°")
|
|
print(f"\nMean R²: {np.mean(torque_segment_scores):.6f}")
|
|
print(f"Median R²: {np.median(torque_segment_scores):.6f}")
|
|
print(f"Std Dev: {np.std(torque_segment_scores):.6f}")
|
|
|
|
# Create histogram of torque R² scores
|
|
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
|
|
ax.hist(torque_segment_scores, bins=30, edgecolor='black', alpha=0.7, color='#9467bd')
|
|
ax.axvline(np.mean(torque_segment_scores), color='red', linestyle='--', linewidth=2, label=f'Mean = {np.mean(torque_segment_scores):.6f}')
|
|
ax.axvline(np.median(torque_segment_scores), color='orange', linestyle='--', linewidth=2, label=f'Median = {np.median(torque_segment_scores):.6f}')
|
|
ax.set_xlabel('R² Score', fontsize=12)
|
|
ax.set_ylabel('Number of Segments', fontsize=12)
|
|
ax.set_title('Distribution of Torque R² Scores Across All Parameter Segments', fontsize=14, fontweight='bold')
|
|
ax.legend(fontsize=11)
|
|
ax.grid(True, alpha=0.3, axis='y')
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
# Get worst torque segment data
|
|
worst_torque_data = magDf.iloc[worst_torque_segment['start']:worst_torque_segment['end']]
|
|
worst_torque_currL = worst_torque_segment['currL']
|
|
worst_torque_currR = worst_torque_segment['currR']
|
|
worst_torque_roll = worst_torque_segment['roll']
|
|
worst_torque_r2 = worst_torque_segment['r2']
|
|
|
|
# Create fine grid for this worst segment
|
|
gap_fine_worst_torque = np.linspace(2, 30, 500)
|
|
X_worst_torque_fine = pd.DataFrame({
|
|
'currL [A]': [worst_torque_currL] * 500,
|
|
'currR [A]': [worst_torque_currR] * 500,
|
|
'rollDeg [deg]': [worst_torque_roll] * 500,
|
|
'invGap': 1/gap_fine_worst_torque
|
|
})
|
|
|
|
# Get predictions
|
|
X_worst_torque_poly = poly_best.transform(X_worst_torque_fine)
|
|
y_worst_torque_pred = model_best.predict(X_worst_torque_poly)
|
|
|
|
# Extract actual data
|
|
gap_worst_torque_actual = worst_torque_data['GapHeight [mm]'].values
|
|
torque_worst_actual = worst_torque_data['YokeTorque.Torque [mNewtonMeter]'].values
|
|
|
|
# Create the plot
|
|
fig, ax = plt.subplots(1, 1, figsize=(12, 7))
|
|
|
|
# Plot polynomial curve for TORQUE
|
|
ax.plot(gap_fine_worst_torque, y_worst_torque_pred[:, 1], '-', linewidth=2.5,
|
|
label=f'Degree {best_degree} Polynomial', color='#ff7f0e', alpha=0.8)
|
|
|
|
# Plot actual data points
|
|
ax.scatter(gap_worst_torque_actual, torque_worst_actual, s=120, marker='o',
|
|
color='#9467bd', edgecolors='black', linewidths=1.5,
|
|
label='Ansys Data', zorder=5)
|
|
|
|
# Formatting
|
|
ax.set_xlabel('Gap Height (mm)', fontsize=13)
|
|
ax.set_ylabel('Torque (mN·m)', fontsize=13)
|
|
ax.set_title(f'Worst Torque Case (R² = {worst_torque_r2:.6f}) - currL={worst_torque_currL:.0f}A, currR={worst_torque_currR:.0f}A, roll={worst_torque_roll:.1f}°',
|
|
fontsize=14, fontweight='bold')
|
|
ax.legend(fontsize=12, loc='best')
|
|
ax.grid(True, alpha=0.3)
|
|
|
|
# Add vertical lines at data points
|
|
for gap_point in gap_worst_torque_actual:
|
|
ax.axvline(gap_point, color='gray', linestyle=':', alpha=0.3, linewidth=0.8)
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print(f"Worst-case torque prediction achieves R² = {worst_torque_r2:.6f}")
|
|
|
|
[12]. Tool: Claude Code
|
|
Prompt:
|
|
Can we now export this model for use in a simulator? I want this to be accessible but also quite
|
|
fast since the simulator will be running at 240Hz and calculating forces and torques for 2 yokes at once.
|
|
|
|
Output:
|
|
print("=" * 60)
|
|
print("MODEL EXPORT SUMMARY")
|
|
print("=" * 60)
|
|
print(f"\nModel Type: Polynomial Regression")
|
|
print(f"Polynomial Degree: {best_degree}")
|
|
print(f"Number of Features: {poly_best.n_features_in_}")
|
|
print(f"Number of Polynomial Terms: {poly_best.n_output_features_}")
|
|
print(f"\nInput Features:")
|
|
print(" - currL [A]")
|
|
print(" - currR [A]")
|
|
print(" - rollDeg [deg]")
|
|
print(" - invGap (1/GapHeight)")
|
|
print(f"\nOutputs:")
|
|
print(" - YokeForce.Force_z [newton]")
|
|
print(" - YokeTorque.Torque [mNewtonMeter]")
|
|
print(f"\nModel Performance:")
|
|
print(f" Force R² Score: {r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]):.6f}")
|
|
print(f" Torque R² Score: {torque_r2_overall:.6f}")
|
|
print(f"\nModel Coefficients Shape:")
|
|
print(f" Force coefficients: {model_best.coef_[0].shape}")
|
|
print(f" Torque coefficients: {model_best.coef_[1].shape}")
|
|
print("=" * 60)
|
|
|
|
# Generate standalone predictor class
|
|
predictor_code = f'''"""
|
|
Magnetic Levitation Force and Torque Predictor
|
|
Generated from polynomial regression model (degree {best_degree})
|
|
|
|
Performance:
|
|
- Force R²: {r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]):.6f}
|
|
- Torque R²: {torque_r2_overall:.6f}
|
|
|
|
Usage:
|
|
predictor = MaglevPredictor()
|
|
force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)
|
|
"""
|
|
|
|
import numpy as np
|
|
from itertools import combinations_with_replacement
|
|
|
|
class MaglevPredictor:
|
|
def __init__(self):
|
|
"""Initialize the magnetic levitation force/torque predictor."""
|
|
self.degree = {best_degree}
|
|
self.n_features = 4 # currL, currR, roll, invGap
|
|
|
|
# Force model coefficients
|
|
self.force_intercept = {model_best.intercept_[0]}
|
|
self.force_coef = np.array({model_best.coef_[0].tolist()})
|
|
|
|
# Torque model coefficients
|
|
self.torque_intercept = {model_best.intercept_[1]}
|
|
self.torque_coef = np.array({model_best.coef_[1].tolist()})
|
|
|
|
def _polynomial_features(self, X):
|
|
"""
|
|
Generate polynomial features up to specified degree.
|
|
Mimics sklearn's PolynomialFeatures with include_bias=True.
|
|
|
|
Args:
|
|
X: numpy array of shape (n_samples, 4) with [currL, currR, roll, invGap]
|
|
|
|
Returns:
|
|
Polynomial features array
|
|
"""
|
|
n_samples = X.shape[0]
|
|
|
|
# Start with bias term (column of ones)
|
|
features = [np.ones(n_samples)]
|
|
|
|
# Add original features
|
|
for i in range(self.n_features):
|
|
features.append(X[:, i])
|
|
|
|
# Add polynomial combinations
|
|
for deg in range(2, self.degree + 1):
|
|
for combo in combinations_with_replacement(range(self.n_features), deg):
|
|
term = np.ones(n_samples)
|
|
for idx in combo:
|
|
term *= X[:, idx]
|
|
features.append(term)
|
|
|
|
return np.column_stack(features)
|
|
|
|
def predict(self, currL, currR, roll, gap_height):
|
|
"""
|
|
Predict force and torque for given operating conditions.
|
|
|
|
Args:
|
|
currL: Left coil current in Amps
|
|
currR: Right coil current in Amps
|
|
roll: Roll angle in degrees
|
|
gap_height: Gap height in mm
|
|
|
|
Returns:
|
|
tuple: (force [N], torque [mN·m])
|
|
"""
|
|
# Compute inverse gap (critical transformation!)
|
|
invGap = 1.0 / gap_height
|
|
|
|
# Create input array
|
|
X = np.array([[currL, currR, roll, invGap]])
|
|
|
|
# Generate polynomial features
|
|
X_poly = self._polynomial_features(X)
|
|
|
|
# Compute predictions
|
|
force = self.force_intercept + np.dot(X_poly, self.force_coef)[0]
|
|
torque = self.torque_intercept + np.dot(X_poly, self.torque_coef)[0]
|
|
|
|
return force, torque
|
|
|
|
def predict_batch(self, currL_array, currR_array, roll_array, gap_height_array):
|
|
"""
|
|
Predict force and torque for multiple operating conditions.
|
|
|
|
Args:
|
|
currL_array: Array of left coil currents [A]
|
|
currR_array: Array of right coil currents [A]
|
|
roll_array: Array of roll angles [deg]
|
|
gap_height_array: Array of gap heights [mm]
|
|
|
|
Returns:
|
|
tuple: (force_array [N], torque_array [mN·m])
|
|
"""
|
|
# Convert to numpy arrays
|
|
currL_array = np.asarray(currL_array)
|
|
currR_array = np.asarray(currR_array)
|
|
roll_array = np.asarray(roll_array)
|
|
gap_height_array = np.asarray(gap_height_array)
|
|
|
|
# Compute inverse gaps
|
|
invGap_array = 1.0 / gap_height_array
|
|
|
|
# Stack into feature matrix
|
|
X = np.column_stack([currL_array, currR_array, roll_array, invGap_array])
|
|
|
|
# Generate polynomial features
|
|
X_poly = self._polynomial_features(X)
|
|
|
|
# Compute predictions
|
|
force_array = self.force_intercept + np.dot(X_poly, self.force_coef)
|
|
torque_array = self.torque_intercept + np.dot(X_poly, self.torque_coef)
|
|
|
|
return force_array, torque_array
|
|
|
|
|
|
if __name__ == "__main__":
|
|
# Example usage
|
|
predictor = MaglevPredictor()
|
|
|
|
# Single prediction
|
|
force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)
|
|
print(f"Single prediction:")
|
|
print(f" Force: {{force:.2f}} N")
|
|
print(f" Torque: {{torque:.2f}} mN·m")
|
|
|
|
# Batch prediction
|
|
currL = np.array([-15, -15, -10])
|
|
currR = np.array([-15, -10, -10])
|
|
roll = np.array([0, 0.5, 1.0])
|
|
gap = np.array([10, 12, 15])
|
|
|
|
forces, torques = predictor.predict_batch(currL, currR, roll, gap)
|
|
print(f"\\nBatch prediction:")
|
|
for i in range(len(forces)):
|
|
print(f" Condition {{i+1}}: Force={{forces[i]:.2f}} N, Torque={{torques[i]:.2f}} mN·m")
|
|
'''
|
|
|
|
# Save to file
|
|
predictor_filename = 'maglev_predictor.py'
|
|
with open(predictor_filename, 'w') as f:
|
|
f.write(predictor_code)
|
|
|
|
print(f"✓ Saved standalone predictor: {predictor_filename}")
|
|
print(f"\nThis module:")
|
|
print(f" - Has NO scikit-learn dependency (only numpy)")
|
|
print(f" - Can be imported directly into your simulator")
|
|
print(f" - Includes both single and batch prediction methods")
|
|
print(f" - Implements polynomial feature generation internally")
|
|
print(f"\nTo use:")
|
|
print(f"```python")
|
|
print(f"from maglev_predictor import MaglevPredictor")
|
|
print(f"predictor = MaglevPredictor()")
|
|
print(f"force, torque = predictor.predict(currL=-15, currR=-15, roll=1.0, gap_height=10.0)")
|
|
print(f"```")
|
|
|
|
# Import the standalone predictor (reload to get latest version)
|
|
import importlib
|
|
import maglev_predictor
|
|
importlib.reload(maglev_predictor)
|
|
from maglev_predictor import MaglevPredictor
|
|
|
|
# Create predictor instance
|
|
standalone_predictor = MaglevPredictor()
|
|
|
|
# Test cases
|
|
test_cases = [
|
|
{'currL': -15, 'currR': -15, 'roll': 0.0, 'gap': 10.0},
|
|
{'currL': -15, 'currR': -15, 'roll': 1.0, 'gap': 10.0},
|
|
{'currL': -15, 'currR': -10, 'roll': 0.5, 'gap': 12.0},
|
|
{'currL': -10, 'currR': -10, 'roll': 2.0, 'gap': 15.0},
|
|
]
|
|
|
|
print("Validation: Comparing Standalone Predictor vs Original Model")
|
|
print("=" * 80)
|
|
|
|
for i, case in enumerate(test_cases):
|
|
# Standalone predictor
|
|
force_standalone, torque_standalone = standalone_predictor.predict(
|
|
case['currL'], case['currR'], case['roll'], case['gap']
|
|
)
|
|
|
|
# Original model
|
|
X_test = pd.DataFrame({
|
|
'currL [A]': [case['currL']],
|
|
'currR [A]': [case['currR']],
|
|
'rollDeg [deg]': [case['roll']],
|
|
'invGap': [1/case['gap']]
|
|
})
|
|
X_test_poly = poly_best.transform(X_test)
|
|
y_test_pred = model_best.predict(X_test_poly)
|
|
force_original = y_test_pred[0, 0]
|
|
torque_original = y_test_pred[0, 1]
|
|
|
|
# Compare
|
|
force_diff = abs(force_standalone - force_original)
|
|
torque_diff = abs(torque_standalone - torque_original)
|
|
|
|
print(f"\nTest Case {i+1}: currL={case['currL']}A, currR={case['currR']}A, " +
|
|
f"roll={case['roll']}°, gap={case['gap']}mm")
|
|
print(f" Force: Standalone={force_standalone:9.4f} N | Original={force_original:9.4f} N | Diff={force_diff:.2e}")
|
|
print(f" Torque: Standalone={torque_standalone:9.4f} mN·m | Original={torque_original:9.4f} mN·m | Diff={torque_diff:.2e}")
|
|
|
|
# Verify match (should be essentially identical, accounting for floating point)
|
|
assert force_diff < 1e-8, f"Force mismatch in test case {i+1}"
|
|
assert torque_diff < 1e-8, f"Torque mismatch in test case {i+1}"
|
|
|
|
print("\n" + "=" * 80)
|
|
print("✓ All tests passed! Standalone predictor matches original model perfectly.")
|
|
print("\nThe standalone predictor is ready for integration into your simulator!")
|
|
|
|
import joblib
|
|
import pickle
|
|
|
|
# Save the polynomial transformer and model
|
|
model_data = {
|
|
'poly_features': poly_best,
|
|
'model': model_best,
|
|
'degree': best_degree,
|
|
'feature_names': ['currL [A]', 'currR [A]', 'rollDeg [deg]', 'invGap'],
|
|
'performance': {
|
|
'force_r2': r2_score(y.iloc[:, 0].values, y_pred_full[:, 0]),
|
|
'torque_r2': torque_r2_overall
|
|
}
|
|
}
|
|
|
|
# Save using joblib (preferred for sklearn models)
|
|
joblib.dump(model_data, 'maglev_model.pkl')
|
|
|
|
print("✓ Saved sklearn model: maglev_model.pkl")
|
|
print(f"\nFile size: {os.path.getsize('maglev_model.pkl') / 1024:.2f} KB")
|
|
print(f"\nContents:")
|
|
print(f" - PolynomialFeatures transformer (degree {best_degree})")
|
|
print(f" - LinearRegression model")
|
|
print(f" - Feature names and metadata")
|
|
print(f" - Performance metrics")
|
|
print(f"\nTo use:")
|
|
print(f"```python")
|
|
print(f"import joblib")
|
|
print(f"import numpy as np")
|
|
print(f"import pandas as pd")
|
|
print(f"")
|
|
print(f"# Load model")
|
|
print(f"data = joblib.load('maglev_model.pkl')")
|
|
print(f"poly = data['poly_features']")
|
|
print(f"model = data['model']")
|
|
print(f"")
|
|
print(f"# Make prediction")
|
|
print(f"X = pd.DataFrame({{'currL [A]': [-15], 'currR [A]': [-15], ")
|
|
print(f" 'rollDeg [deg]': [1.0], 'invGap': [1/10.0]}})")
|
|
print(f"X_poly = poly.transform(X)")
|
|
print(f"force, torque = model.predict(X_poly)[0]")
|
|
print(f"```")
|
|
|
|
import time
|
|
|
|
# Load the sklearn model
|
|
loaded_data = joblib.load('maglev_model.pkl')
|
|
loaded_poly = loaded_data['poly_features']
|
|
loaded_model = loaded_data['model']
|
|
|
|
# Create test data (1000 predictions)
|
|
n_tests = 1000
|
|
test_currL = np.random.uniform(-15, 0, n_tests)
|
|
test_currR = np.random.uniform(-15, 0, n_tests)
|
|
test_roll = np.random.uniform(-5, 5, n_tests)
|
|
test_gap = np.random.uniform(6, 24, n_tests)
|
|
|
|
print("Performance Comparison (1000 predictions)")
|
|
print("=" * 60)
|
|
|
|
# Benchmark 1: Sklearn model (loaded from pickle)
|
|
start = time.perf_counter()
|
|
X_test = pd.DataFrame({
|
|
'currL [A]': test_currL,
|
|
'currR [A]': test_currR,
|
|
'rollDeg [deg]': test_roll,
|
|
'invGap': 1/test_gap
|
|
})
|
|
X_test_poly = loaded_poly.transform(X_test)
|
|
sklearn_results = loaded_model.predict(X_test_poly)
|
|
sklearn_time = time.perf_counter() - start
|
|
|
|
print(f"Sklearn Model (pickled): {sklearn_time*1000:.2f} ms")
|
|
|
|
# Benchmark 2: Standalone predictor
|
|
start = time.perf_counter()
|
|
standalone_results = standalone_predictor.predict_batch(test_currL, test_currR, test_roll, test_gap)
|
|
standalone_time = time.perf_counter() - start
|
|
|
|
print(f"Standalone Predictor: {standalone_time*1000:.2f} ms")
|
|
|
|
# Speedup
|
|
speedup = standalone_time / sklearn_time
|
|
print(f"\nSpeedup: {speedup:.2f}x {'faster' if speedup > 1 else 'slower'} with sklearn model")
|
|
|
|
# Verify results match
|
|
force_diff = np.abs(sklearn_results[:, 0] - standalone_results[0]).max()
|
|
torque_diff = np.abs(sklearn_results[:, 1] - standalone_results[1]).max()
|
|
|
|
print(f"\nResult verification:")
|
|
print(f" Max force difference: {force_diff:.2e} N")
|
|
print(f" Max torque difference: {torque_diff:.2e} mN·m")
|
|
print(f" ✓ Results match!" if (force_diff < 1e-8 and torque_diff < 1e-8) else " ✗ Results differ!")
|
|
|
|
# Per-prediction timing
|
|
print(f"\nPer-prediction time:")
|
|
print(f" Sklearn: {sklearn_time/n_tests*1e6:.2f} μs")
|
|
print(f" Standalone: {standalone_time/n_tests*1e6:.2f} μs")
|
|
|
|
print(f"\n{'='*60}")
|
|
print(f"Recommendation:")
|
|
if speedup > 1.2:
|
|
print(f" Use sklearn model (maglev_model.pkl) - {speedup:.1f}x faster!")
|
|
elif speedup < 0.8:
|
|
print(f" Use standalone predictor - {1/speedup:.1f}x faster!")
|
|
else:
|
|
print(f" Performance is similar - choose based on dependencies:")
|
|
print(f" • Sklearn model: Faster, requires sklearn+joblib")
|
|
print(f" • Standalone: No sklearn dependency, portable")
|
|
|
|
[13]. Tool: Claude Code
|
|
Prompt:
|
|
Now can we use our model to find the equilibrium gab height? ie the gap height at which,
|
|
with currents and roll at 0, our force from a single yoke * 2 = the weight of the pod in newtons?
|
|
|
|
Output:
|
|
# Find equilibrium gap height for 5.8 kg pod using polynomial root finding
|
|
import numpy as np
|
|
from maglev_predictor import MaglevPredictor
|
|
from scipy.optimize import fsolve
|
|
|
|
# Initialize predictor
|
|
predictor = MaglevPredictor()
|
|
|
|
# Target force for 5.8 kg pod (total force = weight)
|
|
# Since we have TWO yokes (front and back), each produces this force
|
|
target_force_per_yoke = 5.8 * 9.81 / 2 # 28.449 N per yoke
|
|
|
|
print("=" * 70)
|
|
print("EQUILIBRIUM GAP HEIGHT FINDER (Analytical Solution)")
|
|
print("=" * 70)
|
|
print(f"Pod mass: 5.8 kg")
|
|
print(f"Total weight: {5.8 * 9.81:.3f} N")
|
|
print(f"Target force per yoke: {target_force_per_yoke:.3f} N")
|
|
print(f"Parameters: currL = 0 A, currR = 0 A, roll = 0°")
|
|
print()
|
|
|
|
# Method 2: Use scipy.optimize.fsolve for verification
|
|
def force_error(gap_height):
|
|
# Handle array input from fsolve (convert to scalar)
|
|
gap_height = float(np.atleast_1d(gap_height)[0])
|
|
force, _ = predictor.predict(currL=0, currR=0, roll=0, gap_height=gap_height)
|
|
return force - target_force_per_yoke
|
|
|
|
# Try multiple initial guesses to find all solutions
|
|
initial_guesses = [8, 10, 15, 20, 25]
|
|
scipy_solutions = []
|
|
|
|
print("using scipy.optimize.fsolve:")
|
|
for guess in initial_guesses:
|
|
solution = fsolve(force_error, guess)[0]
|
|
if 5 <= solution <= 30: # Valid range
|
|
force, torque = predictor.predict(currL=0, currR=0, roll=0, gap_height=solution)
|
|
error = abs(force - target_force_per_yoke)
|
|
if error < 0.01: # Valid solution (within 10 mN)
|
|
scipy_solutions.append((solution, force))
|
|
|
|
# Remove duplicates (solutions within 0.1 mm)
|
|
unique_solutions = []
|
|
for sol, force in scipy_solutions:
|
|
is_duplicate = False
|
|
for unique_sol, _ in unique_solutions:
|
|
if abs(sol - unique_sol) < 0.1:
|
|
is_duplicate = True
|
|
break
|
|
if not is_duplicate:
|
|
unique_solutions.append((sol, force))
|
|
|
|
for gap_val, force in unique_solutions:
|
|
print(f" Gap: {gap_val:.6f} mm → Force: {force:.6f} N")
|
|
print()
|
|
|
|
[14]. Tool: Gemini
|
|
Prompt:
|
|
Can we make a step function that uses our ansys-based regression model to find
|
|
forces on each yoke and apply them to the Bullet simulation? Also, we must
|
|
calculate a reward for any given pod position and orientation, prioritizing
|
|
stable levitation and penalizing sticking to the track or falling off more heavily
|
|
|
|
Output:
|
|
def step(self, action):
|
|
# Check if PyBullet connection is still active (GUI might be closed)
|
|
try:
|
|
p.getConnectionInfo(physicsClientId=self.client)
|
|
except p.error:
|
|
# Connection lost - GUI was closed
|
|
return self._get_obs(), -100.0, True, True, {'error': 'GUI closed'}
|
|
|
|
# Update Coil Currents from PWM Actions
|
|
pwm_front_L = action[0] # yoke +x,+y
|
|
pwm_front_R = action[1] # yoke +x,-y
|
|
pwm_back_L = action[2] # yoke -x,+y
|
|
pwm_back_R = action[3] # yoke -x,-y
|
|
|
|
curr_front_L = self.coil_front_L.update(pwm_front_L, self.dt)
|
|
curr_front_R = self.coil_front_R.update(pwm_front_R, self.dt)
|
|
curr_back_L = self.coil_back_L.update(pwm_back_L, self.dt)
|
|
curr_back_R = self.coil_back_R.update(pwm_back_R, self.dt)
|
|
|
|
# --- 2. Get Current Pod State ---
|
|
pos, orn = p.getBasePositionAndOrientation(self.podId, physicsClientId=self.client)
|
|
lin_vel, ang_vel = p.getBaseVelocity(self.podId, physicsClientId=self.client)
|
|
|
|
# Convert quaternion to rotation matrix
|
|
rot_matrix = np.array(p.getMatrixFromQuaternion(orn)).reshape(3, 3)
|
|
|
|
# --- 3. Calculate Gap Heights at Yoke Positions (for force prediction) ---
|
|
# Calculate world positions of the 4 yokes (ends of U-yokes)
|
|
yoke_gap_heights_dict = {} # Store by label for easy access
|
|
|
|
for i, yoke_idx in enumerate(self.yoke_indices):
|
|
local_pos = self.collision_local_positions[yoke_idx]
|
|
local_vec = np.array(local_pos)
|
|
world_offset = rot_matrix @ local_vec
|
|
world_pos = np.array(pos) + world_offset
|
|
|
|
# Top surface of yoke box (add half-height = 5mm)
|
|
yoke_top_z = world_pos[2] + 0.005
|
|
|
|
# Gap height: track bottom (Z=0) to yoke top (negative Z)
|
|
gap_height = -yoke_top_z # Convert to positive gap in meters
|
|
yoke_gap_heights_dict[self.yoke_labels[i]] = gap_height
|
|
|
|
# Average gap heights for each U-shaped yoke (average left and right ends)
|
|
# Front yoke: average of Front_Left and Front_Right
|
|
# Back yoke: average of Back_Left and Back_Right
|
|
avg_gap_front = (yoke_gap_heights_dict.get('Front_Left', 0.010) +
|
|
yoke_gap_heights_dict.get('Front_Right', 0.010)) / 2
|
|
avg_gap_back = (yoke_gap_heights_dict.get('Back_Left', 0.010) +
|
|
yoke_gap_heights_dict.get('Back_Right', 0.010)) / 2
|
|
|
|
front_left_gap = yoke_gap_heights_dict.get('Front_Left', 0.010)
|
|
front_right_gap = yoke_gap_heights_dict.get('Front_Right', 0.010)
|
|
back_left_gap = yoke_gap_heights_dict.get('Back_Left', 0.010)
|
|
back_right_gap = yoke_gap_heights_dict.get('Back_Right', 0.010)
|
|
|
|
y_distance = 0.1016 # 2 * 0.0508m (left to right distance)
|
|
x_distance = 0.2518 # 2 * 0.1259m (front to back distance)
|
|
|
|
# Roll angle
|
|
# When right side has larger gap, roll is negative
|
|
roll_angle_front = np.arcsin(-(front_right_gap - front_left_gap) / y_distance)
|
|
roll_angle_back = np.arcsin(-(back_right_gap - back_left_gap) / y_distance)
|
|
roll_angle = (roll_angle_front + roll_angle_back) / 2
|
|
|
|
# When back has larger gap, pitch is positive
|
|
pitch_angle_left = np.arcsin((back_left_gap - front_left_gap) / x_distance)
|
|
pitch_angle_right = np.arcsin((back_right_gap - front_right_gap) / x_distance)
|
|
pitch_angle = (pitch_angle_left + pitch_angle_right) / 2
|
|
|
|
# Predict Forces and Torques using Maglev Predictor
|
|
# Gap heights in mm
|
|
gap_front_mm = avg_gap_front * 1000
|
|
gap_back_mm = avg_gap_back * 1000
|
|
|
|
# Roll angle in degrees
|
|
roll_deg = np.degrees(roll_angle)
|
|
|
|
# Predict force and torque for each U-shaped yoke
|
|
# Front yoke
|
|
force_front, torque_front = self.predictor.predict(
|
|
curr_front_L, curr_front_R, roll_deg, gap_front_mm
|
|
)
|
|
|
|
# Back yoke
|
|
force_back, torque_back = self.predictor.predict(
|
|
curr_back_L, curr_back_R, roll_deg, gap_back_mm
|
|
)
|
|
|
|
# --- 5. Apply Forces and Torques to Pod ---
|
|
# Forces are applied at Y=0 (center of U-yoke) at each X position
|
|
# This is where the actual magnetic force acts on the U-shaped yoke
|
|
|
|
# Apply force at front yoke center (X=+0.1259, Y=0)
|
|
front_yoke_center = [0.1259, 0, 0.08585] # From pod.xml yoke positions
|
|
p.applyExternalForce(
|
|
self.podId, -1,
|
|
forceObj=[0, 0, force_front],
|
|
posObj=front_yoke_center,
|
|
flags=p.LINK_FRAME,
|
|
physicsClientId=self.client
|
|
)
|
|
|
|
# Apply force at back yoke center (X=-0.1259, Y=0)
|
|
back_yoke_center = [-0.1259, 0, 0.08585]
|
|
p.applyExternalForce(
|
|
self.podId, -1,
|
|
forceObj=[0, 0, force_back],
|
|
posObj=back_yoke_center,
|
|
flags=p.LINK_FRAME,
|
|
physicsClientId=self.client
|
|
)
|
|
|
|
|
|
|
|
# Apply roll torques
|
|
# Each yoke produces its own torque about X axis
|
|
torque_front_Nm = torque_front / 1000 # Convert from mN·m to N·m
|
|
torque_back_Nm = torque_back / 1000
|
|
|
|
# Apply torques at respective yoke positions
|
|
p.applyExternalTorque(
|
|
self.podId, -1,
|
|
torqueObj=[torque_front_Nm, 0, 0],
|
|
flags=p.LINK_FRAME,
|
|
physicsClientId=self.client
|
|
)
|
|
p.applyExternalTorque(
|
|
self.podId, -1,
|
|
torqueObj=[torque_back_Nm, 0, 0],
|
|
flags=p.LINK_FRAME,
|
|
physicsClientId=self.client
|
|
)
|
|
|
|
# --- 6. Step Simulation ---
|
|
p.stepSimulation(physicsClientId=self.client)
|
|
self.current_step += 1
|
|
|
|
# Check for physical contact with track (bolts touching)
|
|
contact_points = p.getContactPoints(bodyA=self.podId, bodyB=self.trackId, physicsClientId=self.client)
|
|
has_contact = len(contact_points) > 0
|
|
|
|
# --- 7. Get New Observation ---
|
|
obs = self._get_obs()
|
|
|
|
# --- 8. Calculate Reward ---
|
|
# Goal: Hover at target gap (16.5mm), minimize roll/pitch, minimize power
|
|
target_gap = TARGET_GAP # 16.5mm in meters
|
|
avg_gap = (avg_gap_front + avg_gap_back) / 2
|
|
|
|
gap_error = abs(avg_gap - target_gap)
|
|
|
|
# Power dissipation (all 4 coils)
|
|
power = (curr_front_L**2 * self.coil_front_L.R +
|
|
curr_front_R**2 * self.coil_front_R.R +
|
|
curr_back_L**2 * self.coil_back_L.R +
|
|
curr_back_R**2 * self.coil_back_R.R)
|
|
|
|
# --- Improved Reward Function ---
|
|
# Use reward shaping with reasonable scales to enable learning
|
|
|
|
# 1. Gap Error Reward (most important)
|
|
# Use exponential decay for smooth gradient near target
|
|
gap_error_mm = gap_error * 1000 # Convert to mm
|
|
gap_reward = 10.0 * np.exp(-0.5 * (gap_error_mm / 3.0)**2) # Peak at 0mm error, 3mm std dev
|
|
|
|
# 2. Orientation Penalties (smaller scale)
|
|
roll_penalty = abs(np.degrees(roll_angle)) * 0.02
|
|
pitch_penalty = abs(np.degrees(pitch_angle)) * 0.02
|
|
|
|
# 3. Velocity Penalty (discourage rapid oscillations)
|
|
z_velocity = lin_vel[2]
|
|
velocity_penalty = abs(z_velocity) * 0.1
|
|
|
|
# 4. Contact Penalty
|
|
contact_points = p.getContactPoints(bodyA=self.podId, bodyB=self.trackId)
|
|
contact_penalty = len(contact_points) * 0.2
|
|
|
|
# 5. Power Penalty (encourage efficiency, but small weight)
|
|
power_penalty = power * 0.001
|
|
|
|
# Combine rewards (scaled to ~[-5, +1] range per step)
|
|
reward = gap_reward - roll_penalty - pitch_penalty - velocity_penalty - contact_penalty - power_penalty
|
|
|
|
# Check Termination (tighter bounds for safety)
|
|
terminated = False
|
|
truncated = False
|
|
|
|
# Terminate if gap is too small (crash) or too large (lost)
|
|
if avg_gap < 0.003 or avg_gap > 0.035:
|
|
terminated = True
|
|
reward = -10.0 # Failure penalty (scaled down)
|
|
|
|
# Terminate if orientation is too extreme
|
|
if abs(roll_angle) > np.radians(15) or abs(pitch_angle) > np.radians(15):
|
|
terminated = True
|
|
reward = -10.0
|
|
|
|
# Success bonus for stable hovering near target
|
|
if gap_error_mm < 1.0 and abs(np.degrees(roll_angle)) < 2.0 and abs(np.degrees(pitch_angle)) < 2.0:
|
|
reward += 2.0 # Bonus for excellent control
|
|
|
|
info = {
|
|
'curr_front_L': curr_front_L,
|
|
'curr_front_R': curr_front_R,
|
|
'curr_back_L': curr_back_L,
|
|
'curr_back_R': curr_back_R,
|
|
'gap_front_yoke': avg_gap_front,
|
|
'gap_back_yoke': avg_gap_back,
|
|
'roll': roll_angle,
|
|
'force_front': force_front,
|
|
'force_back': force_back,
|
|
'torque_front': torque_front,
|
|
'torque_back': torque_back
|
|
}
|
|
|
|
return obs, reward, terminated, truncated, info
|
|
|
|
[15]. Tool: Gemini
|
|
Prompt:
|
|
Let's also have an observation function that constructs the input vector that the RL agent sees.
|
|
This should consist of the noisy sensor readings and the velocities of each of these sensors.
|
|
They should be normalized in a way that is friendly to a neural network, ie floating point values
|
|
on [-1,1] ideally.
|
|
|
|
Output:
|
|
def _get_obs(self, initial_reset=False):
|
|
"""
|
|
Returns observation: [gaps(4), velocities(4)]
|
|
Uses noisy sensor readings + computed velocities for microcontroller-friendly deployment
|
|
"""
|
|
pos, orn = p.getBasePositionAndOrientation(self.podId, physicsClientId=self.client)
|
|
|
|
# Convert quaternion to rotation matrix
|
|
rot_matrix = np.array(p.getMatrixFromQuaternion(orn)).reshape(3, 3)
|
|
|
|
# Calculate sensor gap heights with noise
|
|
sensor_gap_heights = {}
|
|
|
|
for i, sensor_idx in enumerate(self.sensor_indices):
|
|
local_pos = self.collision_local_positions[sensor_idx]
|
|
local_vec = np.array(local_pos)
|
|
world_offset = rot_matrix @ local_vec
|
|
world_pos = np.array(pos) + world_offset
|
|
|
|
# Top surface of sensor (add half-height = 5mm)
|
|
sensor_top_z = world_pos[2] + 0.005
|
|
|
|
# Gap height: track bottom (Z=0) to sensor top
|
|
gap_height = -sensor_top_z
|
|
|
|
# Add measurement noise
|
|
noisy_gap = gap_height + np.random.normal(0, self.sensor_noise_std)
|
|
sensor_gap_heights[self.sensor_labels[i]] = noisy_gap
|
|
sensor_gap_heights[self.sensor_labels[i]] = gap_height
|
|
|
|
# Pack sensor measurements in consistent order
|
|
# [center_right, center_left, front, back]
|
|
gaps = np.array([
|
|
sensor_gap_heights.get('Center_Right', 0.010),
|
|
sensor_gap_heights.get('Center_Left', 0.010),
|
|
sensor_gap_heights.get('Front', 0.010),
|
|
sensor_gap_heights.get('Back', 0.010)
|
|
], dtype=np.float32)
|
|
|
|
# Compute velocities (d_gap/dt)
|
|
if initial_reset or (self.prev_sensor_gaps is None):
|
|
# First observation - no velocity information yet
|
|
velocities = np.zeros(4, dtype=np.float32)
|
|
else:
|
|
# Compute velocity as finite difference
|
|
velocities = (gaps - self.prev_sensor_gaps) / self.dt
|
|
|
|
# Store for next step
|
|
self.prev_sensor_gaps = gaps.copy()
|
|
|
|
# Normalize observations
|
|
gaps_normalized = (gaps - TARGET_GAP) / self.gap_scale
|
|
velocities_normalized = velocities / self.velocity_scale
|
|
|
|
# Concatenate: [normalized_gaps, normalized_velocities]
|
|
obs = np.concatenate([gaps_normalized, velocities_normalized])
|
|
|
|
return obs
|
|
|
|
[16]. Tool: Gemini
|
|
Prompt:
|
|
I now need a function that will allow the initial height setting to modify
|
|
the position of the bolts in the Bullet sim, mirroring how the rig is set
|
|
up in real life, as detailed in previous chats.
|
|
|
|
Output:
|
|
def _create_modified_urdf(self):
|
|
"""
|
|
Create a modified URDF with bolt positions adjusted based on initial gap height.
|
|
Bolts are at Z = 0.08585 + gap_mm/1000 (relative to pod origin).
|
|
Yokes and sensors remain at Z = 0.08585 (relative to pod origin).
|
|
"""
|
|
import tempfile
|
|
|
|
# Calculate bolt Z position
|
|
bolt_z = 0.08585 + self.initial_gap_mm / 1000.0
|
|
|
|
# Read the original URDF template
|
|
urdf_template_path = os.path.join(os.path.dirname(__file__), "pod.xml")
|
|
with open(urdf_template_path, 'r') as f:
|
|
urdf_content = f.read()
|
|
|
|
# Replace the bolt Z positions (originally at 0.09585)
|
|
# There are 4 bolts at different X,Y positions but same Z
|
|
urdf_modified = urdf_content.replace(
|
|
'xyz="0.285 0.03 0.09585"',
|
|
f'xyz="0.285 0.03 {bolt_z:.6f}"'
|
|
).replace(
|
|
'xyz="0.285 -0.03 0.09585"',
|
|
f'xyz="0.285 -0.03 {bolt_z:.6f}"'
|
|
).replace(
|
|
'xyz="-0.285 0.03 0.09585"',
|
|
f'xyz="-0.285 0.03 {bolt_z:.6f}"'
|
|
).replace(
|
|
'xyz="-0.285 -0.03 0.09585"',
|
|
f'xyz="-0.285 -0.03 {bolt_z:.6f}"'
|
|
)
|
|
|
|
# Write to a temporary file
|
|
with tempfile.NamedTemporaryFile(mode='w', suffix='.urdf', delete=False) as f:
|
|
f.write(urdf_modified)
|
|
temp_urdf_path = f.name
|
|
|
|
return temp_urdf_path
|
|
|
|
[17]. Tool: Gemini
|
|
Prompt:
|
|
Can you generate a python class that simulates the coil inductance and resistance,
|
|
taking in a PWM duty cycle and applying a voltage based on that duty cycle to get
|
|
an output current?
|
|
|
|
Output:
|
|
class MagLevCoil:
|
|
def __init__(self, r_resistance, l_inductance, source_voltage, maxCurrent):
|
|
self.R = r_resistance
|
|
self.L = l_inductance
|
|
self.current = 0.0
|
|
self.Vs = source_voltage
|
|
self.Imax = maxCurrent
|
|
|
|
def update(self, pwm_duty_cycle, dt):
|
|
"""
|
|
Simulates the coil circuit and force generation.
|
|
pwm_duty_cycle: -1.0 to 1.0
|
|
"""
|
|
# Simple First-order RL circuit approximation
|
|
# V_in = Duty * V_source
|
|
voltage = pwm_duty_cycle * self.Vs
|
|
|
|
# di/dt = (V - I*R) / L
|
|
di = (voltage - self.current * self.R) / self.L
|
|
self.current += di * dt
|
|
self.current = min(max(-self.Imax, self.current), self.Imax)
|
|
|
|
return self.current
|
|
|
|
[18]. Tool: Claude Code
|
|
Prompt:
|
|
Can you generate a testing script that allows me to input a starting height
|
|
and dummy PWM values to the levitation environment so I can do a sanity check
|
|
and see if the environment behaves correctly?
|
|
|
|
Output:
|
|
"""
|
|
Test script for LevPodEnv
|
|
Runs a simple episode with constant actions to verify the environment works
|
|
"""
|
|
|
|
from lev_pod_env import LevPodEnv
|
|
import numpy as np
|
|
import time
|
|
|
|
# Create environment with GUI for visualization
|
|
env = LevPodEnv(use_gui=True, initial_gap_mm=15)
|
|
|
|
print("=" * 60)
|
|
print("Testing LevPodEnv")
|
|
print("=" * 60)
|
|
print(f"Action space: {env.action_space}")
|
|
print(f" 4 PWM duty cycles: [front_L, front_R, back_L, back_R]")
|
|
print(f"Observation space: {env.observation_space}")
|
|
print(f" 8 values: [gaps(4), velocities(4)]")
|
|
print("=" * 60)
|
|
|
|
# Reset environment
|
|
obs, info = env.reset()
|
|
print(f"\nInitial observation:")
|
|
print(f" Gaps: CR={obs[0]*1000:.2f}mm, CL={obs[1]*1000:.2f}mm, F={obs[2]*1000:.2f}mm, B={obs[3]*1000:.2f}mm")
|
|
print(f" Velocities: CR={obs[4]*1000:.2f}mm/s, CL={obs[5]*1000:.2f}mm/s, F={obs[6]*1000:.2f}mm/s, B={obs[7]*1000:.2f}mm/s")
|
|
print(f" Average gap: {np.mean(obs[:4])*1000:.2f} mm")
|
|
|
|
# Run a few steps with constant action to test force application
|
|
print("\nRunning test episode...")
|
|
for step in range(500):
|
|
# Apply constant moderate PWM to all 4 coils
|
|
# 50% PWM should generate current that produces upward force
|
|
action = np.array([0,0,0,0], dtype=np.float32)
|
|
|
|
obs, reward, terminated, truncated, info = env.step(action)
|
|
|
|
if step % 5 == 0:
|
|
print(f"\nStep {step} (t={step/240:.2f}s):")
|
|
print(f" Sensor gaps: CR={obs[0]*1000:.2f}mm, CL={obs[1]*1000:.2f}mm, " +
|
|
f"F={obs[2]*1000:.2f}mm, B={obs[3]*1000:.2f}mm")
|
|
print(f" Velocities: CR={obs[4]*1000:.2f}mm/s, CL={obs[5]*1000:.2f}mm/s, " +
|
|
f"F={obs[6]*1000:.2f}mm/s, B={obs[7]*1000:.2f}mm/s")
|
|
print(f" Yoke gaps: front={info['gap_front_yoke']*1000:.2f}mm, back={info['gap_back_yoke']*1000:.2f}mm")
|
|
print(f" Roll: {np.degrees(info['roll']):.2f}°")
|
|
print(f" Currents: FL={info['curr_front_L']:.2f}A, FR={info['curr_front_R']:.2f}A, " +
|
|
f"BL={info['curr_back_L']:.2f}A, BR={info['curr_back_R']:.2f}A")
|
|
print(f" Forces: front={info['force_front']:.2f}N, back={info['force_back']:.2f}N")
|
|
print(f" Torques: front={info['torque_front']:.2f}mN·m, back={info['torque_back']:.2f}mN·m")
|
|
print(f" Reward: {reward:.2f}")
|
|
|
|
if terminated or truncated:
|
|
print(f"\nEpisode terminated at step {step}")
|
|
break
|
|
|
|
time.sleep(0.01)
|
|
|
|
print("\n" + "=" * 60)
|
|
print("Test complete!")
|
|
print("=" * 60)
|
|
|
|
env.close()
|
|
|
|
[19]. Tool: Claude Code
|
|
Prompt:
|
|
Can we intialize the existing model structure in the rest of the file
|
|
around our use case, with the environment being the lev_pod_env? Additionally,
|
|
I want to be able to see what the theoretical maximum reward is, so I can see how
|
|
well the agent is performing over time without having to use the GUI simulation
|
|
environment.
|
|
|
|
Output:
|
|
environ = LevPodEnv(use_gui=False, initial_gap_mm=14, max_steps=500) # Start below target
|
|
model = ActorCriticNetwork(environ.observation_space.shape[0], environ.action_space.shape[0]).to(device)
|
|
train_data, reward, gap_error = rollout(model, environ)
|
|
print("Reward:", reward)
|
|
print("Episode length:", len(train_data[0]))
|
|
print("Reward per step:", reward / len(train_data[0]))
|
|
print(f"Average gap error: {gap_error:.3f} mm")
|
|
|
|
# Calculate theoretical maximum reward for reference (NEW REWARD STRUCTURE)
|
|
max_steps = 500
|
|
# Best case: gap_reward=1.0, no penalties, success_bonus=2.0
|
|
theoretical_max_per_step = 3.0 # 1.0 (gap) + 2.0 (success bonus) + 0 (no penalties)
|
|
# Realistic good case: gap_reward~0.9, small penalties
|
|
realistic_good_per_step = 0.8
|
|
theoretical_max_total = theoretical_max_per_step * max_steps
|
|
realistic_good_total = realistic_good_per_step * max_steps
|
|
|
|
print(f"\nTheoretical Maximum:")
|
|
print(f" Best per step: {theoretical_max_per_step:.2f}")
|
|
print(f" Realistic good per step: {realistic_good_per_step:.2f}")
|
|
print(f" Best total (500 steps): {theoretical_max_total:.0f}")
|
|
print(f" Current % of realistic: {(reward/realistic_good_total)*100:.1f}%")
|
|
|
|
[20]. Tool: Claude Code
|
|
Prompt:
|
|
Can we improve our reward function, and see how it drops off as the gap error increases?
|
|
It seems the pod is getting too comfortable sitting at the track or just falling away and
|
|
not learning to levitate at a very specific gap height.
|
|
|
|
Output:
|
|
# Visualize the new reward function
|
|
gap_errors_mm = np.linspace(0, 20, 100)
|
|
gap_rewards = np.exp(-0.5 * (gap_errors_mm / 3.0)**2)
|
|
|
|
plt.figure(figsize=(10, 6))
|
|
plt.plot(gap_errors_mm, gap_rewards, linewidth=2)
|
|
plt.axvline(x=1.0, color='g', linestyle='--', label='Success bonus threshold (1mm)')
|
|
plt.axvline(x=3.0, color='orange', linestyle='--', label='1 std dev (3mm)')
|
|
plt.xlabel('Gap Error (mm)', fontsize=12)
|
|
plt.ylabel('Gap Reward Component', fontsize=12)
|
|
plt.title('New Reward Function: Smooth Exponential Decay', fontsize=14, fontweight='bold')
|
|
plt.grid(True, alpha=0.3)
|
|
plt.legend()
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
print("Reward examples at different gap errors:")
|
|
for err in [0, 1, 3, 5, 10, 15]:
|
|
reward = np.exp(-0.5 * (err / 3.0)**2)
|
|
print(f" {err:2d}mm error → {reward:.4f} reward")
|
|
print("\nCompare to old reward: (error*100)² would be:")
|
|
for err in [1, 3, 5, 10]:
|
|
old_penalty = (err * 100)**2
|
|
print(f" {err:2d}mm error → -{old_penalty:,} penalty")
|
|
|
|
[21]. Tool: Claude Code
|
|
Prompt:
|
|
I want to store log files that show the accumulated reward over time
|
|
for each trial, as well as a mean gap error chart that should show whether
|
|
the model is learning to control at 0 gap error.
|
|
|
|
Output:
|
|
import os
|
|
from datetime import datetime
|
|
|
|
# Create RL_Trials folder if it doesn't exist
|
|
os.makedirs('RL_Trials', exist_ok=True)
|
|
|
|
# Create timestamped log file
|
|
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
|
log_file_path = f'RL_Trials/training_log_{timestamp}.txt'
|
|
plot_file_path = f'RL_Trials/gap_error_plot_{timestamp}.png'
|
|
|
|
# Define training params
|
|
num_episodes = 2000 # Increased for more learning time
|
|
print_freq = 20 # Print less frequently
|
|
gui_freq = 100 # Show GUI every 100 episodes
|
|
|
|
# Create PPO trainer with improved hyperparameters
|
|
ppo = PPOTrainer(
|
|
model,
|
|
policy_lr=5e-4, # Higher learning rate
|
|
value_lr=1e-3, # Even higher for value function
|
|
target_kl_div=0.02, # Allow more policy updates
|
|
max_policy_train_iters=40,
|
|
value_train_iters=40,
|
|
entropy_coef=0.02 # More exploration
|
|
)
|
|
|
|
# Open log file
|
|
log_file = open(log_file_path, 'w')
|
|
log_file.write(f"Training Started: {timestamp}\n")
|
|
log_file.write(f"Number of Episodes: {num_episodes}\n")
|
|
log_file.write(f"Print Frequency: {print_freq}\n")
|
|
log_file.write(f"Target Gap Height: {16.491741} mm\n")
|
|
log_file.write(f"Network: 256 hidden units with LayerNorm\n")
|
|
log_file.write(f"Policy LR: 5e-4, Value LR: 1e-3, Entropy: 0.02\n")
|
|
log_file.write("="*70 + "\n\n")
|
|
log_file.flush()
|
|
|
|
print(f"Logging to: {log_file_path}")
|
|
print(f"Plot will be saved to: {plot_file_path}")
|
|
|
|
# Create and save gap error plot
|
|
plt.figure(figsize=(12, 6))
|
|
plt.plot(ep_gap_errors, alpha=0.3, label='Per Episode')
|
|
|
|
# Calculate moving average (window size = print_freq)
|
|
window_size = print_freq
|
|
if len(ep_gap_errors) >= window_size:
|
|
moving_avg = np.convolve(ep_gap_errors, np.ones(window_size)/window_size, mode='valid')
|
|
plt.plot(range(window_size-1, len(ep_gap_errors)), moving_avg, linewidth=2, label=f'{window_size}-Episode Moving Average')
|
|
|
|
plt.xlabel('Episode', fontsize=12)
|
|
plt.ylabel('Average Gap Height Error (mm)', fontsize=12)
|
|
plt.title('Gap Height Error Over Training, n_steps=500', fontsize=14, fontweight='bold')
|
|
plt.legend(fontsize=11)
|
|
plt.grid(True, alpha=0.3)
|
|
plt.tight_layout()
|
|
plt.savefig(plot_file_path, dpi=150, bbox_inches='tight')
|
|
plt.show()
|
|
|
|
print(f"\n✓ Training log saved to: {log_file_path}")
|
|
print(f"✓ Gap error plot saved to: {plot_file_path}") |