From 202a4d7df4071f5a547f44e96f18efba51b9e14f Mon Sep 17 00:00:00 2001 From: DO-Ui <68296986+DO-Ui@users.noreply.github.com> Date: Thu, 19 Feb 2026 13:24:08 -0500 Subject: [PATCH 1/3] Numba Optimizations --- cvtModel/README.md | 4 + cvtModel/setup.py | 17 +- .../models/pulley/primary_pulley_flyweight.py | 28 +- .../models/pulley/pulley_interface.py | 28 +- .../secondary_pulley_torque_reactive.py | 53 ++- .../src/cvt_simulator/models/slip_model.py | 57 +-- .../src/cvt_simulator/simulation_runner.py | 325 +++++++++++++++--- .../src/cvt_simulator/utils/numba_kernels.py | 171 +++++++++ 8 files changed, 513 insertions(+), 170 deletions(-) create mode 100644 cvtModel/src/cvt_simulator/utils/numba_kernels.py diff --git a/cvtModel/README.md b/cvtModel/README.md index bcc26ce9..9eff9ccc 100644 --- a/cvtModel/README.md +++ b/cvtModel/README.md @@ -9,6 +9,10 @@ venv\Scripts\activate # or 'source venv/bin/activate' on Mac ```bash pip install -e .[dev] ``` +For optional runtime acceleration kernels: +```bash +pip install -e .[dev,numba] +``` 3. Run the tests ```bash coverage run -m unittest discover -s test/simulations -s test/utils diff --git a/cvtModel/setup.py b/cvtModel/setup.py index 5d4c8809..83c9baeb 100644 --- a/cvtModel/setup.py +++ b/cvtModel/setup.py @@ -12,10 +12,15 @@ "pandas", "pydantic" ], - extras_require={"dev": [ - "black", - "flake8", - "pytest", - "coverage" - ]}, + extras_require={ + "dev": [ + "black", + "flake8", + "pytest", + "coverage" + ], + "numba": [ + "numba" + ] + }, ) diff --git a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py index 15e3e733..03c90e08 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py +++ b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py @@ -18,6 +18,10 @@ INITIAL_FLYWEIGHT_RADIUS, ) from cvt_simulator.utils.system_state import SystemState +from cvt_simulator.utils.numba_kernels import ( + primary_flyweight_force_kernel, + max_torque_primary_kernel, +) # TODO: Remove this code @@ -152,14 +156,7 @@ def calculate_max_torque( wrap_angle = self._get_wrap_angle(state.shift_distance) radius = self._get_radius(state.shift_distance) - # Capstan equation with V-belt friction enhancement - exp_term = math.exp(self.μ * wrap_angle) - capstan_term = (exp_term - 1) / (exp_term + 1) - radial_force_term = total_radial * radius / np.sin(wrap_angle / 2) - - max_torque = capstan_term * radial_force_term - - return max(0.0, max_torque) # Ensure non-negative + return max_torque_primary_kernel(self.μ, wrap_angle, total_radial, radius) # Private helper methods for force calculations @@ -177,28 +174,19 @@ def _calculate_flyweight_force( shift_distance ) - # Centrifugal force on flyweight: F = m * ω² * r - centrifugal_force = tm.centrifugal_force( + angle, centrifugal_force, angle_multiplier, net = primary_flyweight_force_kernel( self.flyweight_mass, angular_velocity, flyweight_radius, + self.ramp.slope(shift_distance), ) - # Ramp angle at current position - # Ramp is default negative slope, so negate for angle - # If a positive slope ramp is passed, it will generate a force against shifting - # which is expected for such a stupid ramp design. - angle = np.arctan(-self.ramp.slope(shift_distance)) - - # Convert centrifugal force to axial clamping force through ramp angle - net = centrifugal_force * np.tan(angle) - return flyweightForceBreakdown( radius=flyweight_radius, angular_velocity=angular_velocity, angle=angle, centrifugal_force=centrifugal_force, - angle_multiplier=np.tan(angle), + angle_multiplier=angle_multiplier, net=net, ) diff --git a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py index 6b0b0ca5..066ff584 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py +++ b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py @@ -35,6 +35,7 @@ ) from cvt_simulator.models.dataTypes import PulleyState, PulleyForces, PulleyBreakdowns from cvt_simulator.utils.theoretical_models import TheoreticalModels as tm +from cvt_simulator.utils.numba_kernels import radial_force_kernel def get_kwarg(kwargs: dict[str, Any], key: str, default: Any = None) -> Any: @@ -172,11 +173,6 @@ def calculate_radial_force( # radius = self._get_radius(state.shift_distance) # angular_velocity = self._get_angular_velocity(state) - # Radial force from pulley clamping (through V-belt wedging) - radial_from_clamping = ( - 2 * (clamping_force * np.tan(SHEAVE_ANGLE / 2)) / wrap_angle - ) - # To future Kai: # Trust me! If you're using your current assumption # (Which the belt isn't slipping from the secondary, i.e. belt speed = sec speed) @@ -187,18 +183,16 @@ def calculate_radial_force( wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS sec_angular_velocity = state.car_velocity * wheel_to_sec_ratio - radial_from_centrifugal = ( - sec_angular_velocity**2 - * sec_radius**2 - * BELT_CROSS_SECTIONAL_AREA - * RUBBER_DENSITY - ) - - # Total radial force (determines friction capacity) - total_radial = ( - 2 - * np.sin(wrap_angle / 2) - * (radial_from_clamping + radial_from_centrifugal) + radial_from_clamping, radial_from_centrifugal, total_radial = ( + radial_force_kernel( + clamping_force, + SHEAVE_ANGLE, + wrap_angle, + sec_angular_velocity, + sec_radius, + BELT_CROSS_SECTIONAL_AREA, + RUBBER_DENSITY, + ) ) return radial_from_clamping, radial_from_centrifugal, total_radial diff --git a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py index 9c4faab8..333070c3 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py +++ b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py @@ -16,6 +16,10 @@ ) from cvt_simulator.models.ramps import LinearSegment, PiecewiseRamp from cvt_simulator.utils.system_state import SystemState +from cvt_simulator.utils.numba_kernels import ( + secondary_helix_force_kernel, + max_torque_secondary_kernel, +) # TODO: Remove this code @@ -158,36 +162,25 @@ def calculate_max_torque( helix_breakdown = self._calculate_helix_force(0, shift_distance) spring_tors_torque = helix_breakdown.net - # Convert to radial force contribution - spring_force_term = (spring_comp_force + spring_tors_torque) * np.tan( - SHEAVE_ANGLE / 2 - ) - # Centrifugal force contribution, used built in calc with 0 clamp (since we only need centrifugal) _, radial_from_centrifugal, _ = self.calculate_radial_force(state, 0) - centrifugal_force = radial_from_centrifugal * wrap_angle / 2 - - # Capstan term - exp_term = np.exp(self.μ * wrap_angle) - capstan_term = (wrap_angle / (4 * radius)) * (exp_term + 1) / (exp_term - 1) # Torque transmission term (feedback loop) cvt_ratio = tm.current_cvt_ratio(shift_distance) helix_angle = helix_breakdown.angle - transmission_term = ( - 2 - * cvt_ratio - * (HELIX_RADIUS * np.tan(helix_angle)) - * np.tan(SHEAVE_ANGLE / 2) + return max_torque_secondary_kernel( + spring_comp_force, + spring_tors_torque, + SHEAVE_ANGLE, + radial_from_centrifugal, + wrap_angle, + radius, + self.μ, + cvt_ratio, + HELIX_RADIUS, + helix_angle, ) - # Solve for max torque (equilibrium of torque feedback loop) - numerator = centrifugal_force + spring_force_term - denominator = capstan_term - transmission_term - max_torque = numerator / denominator - - return max(0.0, max_torque) # Ensure non-negative - # Private helper methods for force calculations def _calculate_helix_force( @@ -208,16 +201,12 @@ def _calculate_helix_force( # Effective radius at current shift position secondary_radius = tm.outer_sec_radius(shift_distance) - BELT_HEIGHT / 2 - # Helix angle at current position - # Negative slope because helix ramps down with increasing shift - helix_angle = np.arctan(-self.ramp.slope(shift_distance)) - - # Convert torque to axial force through helix geometry - # F = (τ + τ_spring) / (2 * tan(α) * r) - angle_multiplier = 2 * np.tan(helix_angle) * secondary_radius - - # Net - net = (torque + spring_torque_breakdown.net) / angle_multiplier + helix_angle, angle_multiplier, net = secondary_helix_force_kernel( + torque, + spring_torque_breakdown.net, + secondary_radius, + self.ramp.slope(shift_distance), + ) return HelixForceBreakdown( feedbackTorque=torque, diff --git a/cvtModel/src/cvt_simulator/models/slip_model.py b/cvtModel/src/cvt_simulator/models/slip_model.py index 45062271..7d72df82 100644 --- a/cvtModel/src/cvt_simulator/models/slip_model.py +++ b/cvtModel/src/cvt_simulator/models/slip_model.py @@ -16,6 +16,11 @@ from cvt_simulator.constants.constants import ( RUBBER_ALUMINUM_STATIC_FRICTION, ) +from cvt_simulator.utils.numba_kernels import ( + slip_relative_speed_kernel, + slip_coupling_torque_kernel, + torque_demand_kernel, +) class SlipModel: @@ -57,29 +62,18 @@ def get_breakdown(self, state: SystemState) -> SlipBreakdown: torque_demand = self.get_torque_demand(state) wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS - relative_speed = self._relative_speed( + relative_speed = slip_relative_speed_kernel( state.engine_angular_velocity, state.car_velocity * wheel_to_sec_ratio, tm.current_cvt_ratio(state.shift_distance), ) - - # 1) Smooth Coulomb-like torque based on slip speed - # For large |relative_speed|, tanh -> ±1, so |coupling_torque| -> t_max_capacity - # For small |relative_speed|, torque ~ (t_max_capacity / slip_speed_smoothing) * relative_speed (viscous-ish) - coulomb_torque = t_max_capacity * np.tanh( - relative_speed / self.slip_speed_smoothing + coupling_torque, _ = slip_coupling_torque_kernel( + relative_speed, + torque_demand, + t_max_capacity, + self.slip_speed_smoothing, ) - # 2) Optionally respect torque_demand near zero slip by blending - # When relative_speed is small, use torque_demand (clamped); - # as slip grows, fade to the Coulomb model. - v_blend = self.slip_speed_smoothing # you can use same scale or a separate one - alpha = np.clip(abs(relative_speed) / v_blend, 0.0, 1.0) - - torque_demand_clamped = np.clip(torque_demand, -t_max_capacity, t_max_capacity) - - coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque - # 3) Define is_slipping for diagnostics (no effect on dynamics) is_slipping = abs(relative_speed) > self.slip_speed_threshold @@ -111,21 +105,16 @@ def get_torque_demand(self, state: SystemState): * GEARBOX_RATIO ) - # The secret butter - eng_term = engine_torque * wheel_inertia - load_term = ENGINE_INERTIA * load_torque * engine_to_wheel_ratio - shift_term = ( - ENGINE_INERTIA - * wheel_inertia - * wheel_angular_velocity - * engine_to_wheel_ratio_rate_of_change + coupling_torque = torque_demand_kernel( + engine_torque, + load_torque, + wheel_inertia, + wheel_angular_velocity, + engine_to_wheel_ratio, + engine_to_wheel_ratio_rate_of_change, + ENGINE_INERTIA, ) - numerator = eng_term + load_term - shift_term - denominator = wheel_inertia + ENGINE_INERTIA * engine_to_wheel_ratio**2 - - coupling_torque = numerator / denominator - return coupling_torque def get_wheel_speed(self, car_velocity: float): @@ -152,14 +141,6 @@ def calculate_t_max(self, state: SystemState) -> tuple[float, float]: secondary_t_max = max(0, secondary_t_max) return primary_t_max, secondary_t_max - def _relative_speed( - self, - primary_angular_velocity: float, - secondary_angular_velocity: float, - cvt_ratio: float, - ) -> float: - return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) - # def _is_slipping( # self, # primary_angular_velocity: float, diff --git a/cvtModel/src/cvt_simulator/simulation_runner.py b/cvtModel/src/cvt_simulator/simulation_runner.py index ea284c42..ac202cf0 100644 --- a/cvtModel/src/cvt_simulator/simulation_runner.py +++ b/cvtModel/src/cvt_simulator/simulation_runner.py @@ -1,4 +1,5 @@ import sys +import math import numpy as np from typing import Callable, Optional from scipy.integrate import solve_ivp @@ -9,14 +10,16 @@ GEARBOX_RATIO, WHEEL_RADIUS, MAX_SHIFT, + ENGINE_INERTIA, + DRIVELINE_INERTIA, ) from cvt_simulator.utils.conversions import rpm_to_rad_s from cvt_simulator.utils.theoretical_models import TheoreticalModels as tm -from cvt_simulator.utils.simulation_constraints import ( - car_velocity_constraint_event, - get_shift_steady_event, - get_back_shift_event, - shift_constraint_event, +from cvt_simulator.constants.engine_specs import safe_torque_curve +from cvt_simulator.utils.numba_kernels import ( + slip_relative_speed_kernel, + slip_coupling_torque_kernel, + torque_demand_kernel, ) @@ -52,6 +55,63 @@ def __init__( self.progress_callback = progress_callback self._last_callback_percent = -1.0 + # Cached model references for faster ODE RHS evaluation + self._slip_model = system_model.slip_model + self._cvt_shift_model = system_model.cvt_shift_model + self._primary_pulley = self._slip_model.primary_pulley + self._secondary_pulley = self._slip_model.secondary_pulley + self._engine_model = self._slip_model.engine_model + self._load_model = self._slip_model.load_model + + self._wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS + self._wheel_inertia = DRIVELINE_INERTIA + self._slip_model.car_mass * ( + WHEEL_RADIUS**2 + ) + self._cvt_moving_mass = self._cvt_shift_model.cvt_moving_mass + self._slip_speed_smoothing = self._slip_model.slip_speed_smoothing + + self._incline_force = ( + self._load_model.car_mass + * self._load_model.g + * math.sin(self._load_model.incline_angle) + ) + self._drag_force_coeff = ( + 0.5 + * self._load_model.air_density + * self._load_model.frontal_area + * self._load_model.drag_coefficient + ) + + # Precompute CVT ratio and derivative lookup tables to avoid per-step root finding + self._ratio_lut_shift = np.linspace(0.0, MAX_SHIFT, 1024) + self._ratio_lut = np.array( + [tm.current_cvt_ratio(d) for d in self._ratio_lut_shift], dtype=float + ) + self._ratio_rate_per_vel_lut = np.array( + [ + tm.current_cvt_ratio_rate_of_change(d, 1.0) + for d in self._ratio_lut_shift + ], + dtype=float, + ) + + # Precompute torque lookup table to avoid repeated scipy interp1d overhead + self._torque_lut_omega = np.linspace(0.0, rpm_to_rad_s(6000), 2048) + self._torque_lut = np.asarray(safe_torque_curve(self._torque_lut_omega), dtype=float) + + # Reused mutable state object to minimize allocation churn + self._scratch_state = SystemState() + + # Secondary ramp can have a slightly smaller usable max than MAX_SHIFT. + # Clamp internal pulley calculations to avoid inverse-height edge errors. + secondary_last_segment = self._secondary_pulley.ramp.segments[-1] + secondary_max_shift = abs( + self._secondary_pulley.ramp.height(secondary_last_segment.x_end) + ) + self._pulley_calc_shift_max = min(MAX_SHIFT, secondary_max_shift) + + + def run_simulation(self) -> SimulationResult: """Run the simulation and return results.""" cvt_system_ode = self._get_ode_function() @@ -67,9 +127,9 @@ def run_simulation(self) -> SimulationResult: # Phase 1: Normal shifting until full shift is reached events_phase1 = [ - get_shift_steady_event(self.system_model), - car_velocity_constraint_event, - shift_constraint_event, + self._get_shift_steady_event(), + self._get_car_velocity_constraint_event(), + self._get_shift_constraint_event(), ] time_eval_phase1 = time_eval[time_eval >= current_time] @@ -102,8 +162,8 @@ def run_simulation(self) -> SimulationResult: # At full shift, check for back-shift event events_phase2 = [ - get_back_shift_event(self.system_model), - car_velocity_constraint_event, + self._get_back_shift_event(), + self._get_car_velocity_constraint_event(), ] solution_phase2 = self._solve( @@ -129,9 +189,9 @@ def run_simulation(self) -> SimulationResult: break events_phase3 = [ - get_shift_steady_event(self.system_model), - car_velocity_constraint_event, - shift_constraint_event, + self._get_shift_steady_event(), + self._get_car_velocity_constraint_event(), + self._get_shift_constraint_event(), ] solution_phase3 = self._solve( @@ -213,73 +273,224 @@ def _print_progress(self, t): self._last_callback_percent = rounded_percent self.progress_callback(progress_percent) - def _evaluate_cvt_system(self, t: float, y: list[float]): - """Evaluate system dynamics (phase 1: not at full shift).""" - state = SystemState.from_array(y) - self._print_progress(t) + def _get_shift_constraint_event(self): + def shift_constraint_event(t, y): + shift_velocity = y[2] + shift_distance = y[3] + + if shift_distance < 0: + y[3] = 0.0 + y[2] = max(0.0, shift_velocity) + elif shift_distance > MAX_SHIFT: + y[3] = MAX_SHIFT + y[2] = min(0.0, shift_velocity) + + return 1.0 + + return shift_constraint_event + + def _get_car_velocity_constraint_event(self): + def car_velocity_event(t, y): + return y[0] + + car_velocity_event.terminal = True + car_velocity_event.direction = -1 + return car_velocity_event + + def _get_shift_steady_event(self): + def shift_steady_event(t, y): + tol = 1e-5 + if y[3] < MAX_SHIFT - tol: + return -tol + + shift_velocity = y[2] + shift_distance = y[3] + if shift_distance < 0: + y[3] = 0.0 + y[2] = max(0.0, shift_velocity) + elif shift_distance > MAX_SHIFT: + y[3] = MAX_SHIFT + y[2] = min(0.0, shift_velocity) + + _, _, shift_acceleration = self._compute_dynamics(y, full_shift=False) + return shift_acceleration + + shift_steady_event.terminal = True + shift_steady_event.direction = 1 + return shift_steady_event + + def _get_back_shift_event(self): + def back_shift_event(t, y): + if y[3] < MAX_SHIFT - 1e-5: + return 1.0 + + _, _, shift_acceleration = self._compute_dynamics(y, full_shift=False) + return shift_acceleration + 5.0 + + back_shift_event.terminal = True + back_shift_event.direction = -1 + return back_shift_event + + def _clamp_shift_state(self, y: list[float], full_shift: bool) -> tuple[float, float]: + if full_shift: + y[3] = MAX_SHIFT + y[2] = 0.0 + return MAX_SHIFT, 0.0 + + shift_velocity = y[2] + shift_distance = y[3] - # TODO: Remove this (should be handled by constraints) - shift_velocity = state.shift_velocity - shift_distance = state.shift_distance if shift_distance <= 0: - state.shift_distance = 0 - state.shift_velocity = max(0, shift_velocity) - + y[3] = 0.0 + y[2] = max(0.0, shift_velocity) elif shift_distance > MAX_SHIFT: - state.shift_distance = MAX_SHIFT - state.shift_velocity = min(0, shift_velocity) + y[3] = MAX_SHIFT + y[2] = min(0.0, shift_velocity) + + return y[3], y[2] - constrained_y = state.to_array() - for i in range(len(y)): - y[i] = constrained_y[i] + def _lookup_cvt_ratio(self, shift_distance: float) -> float: + return float(np.interp(shift_distance, self._ratio_lut_shift, self._ratio_lut)) - # Get system breakdown (this calculates everything in correct order) - system_breakdown = self.system_model.get_breakdown(state) + def _lookup_cvt_ratio_rate(self, shift_distance: float, shift_velocity: float) -> float: + rate_per_velocity = float( + np.interp(shift_distance, self._ratio_lut_shift, self._ratio_rate_per_vel_lut) + ) + return rate_per_velocity * shift_velocity + + def _lookup_engine_torque(self, engine_angular_velocity: float) -> float: + if engine_angular_velocity <= self._torque_lut_omega[0]: + return float(self._torque_lut[0]) + if engine_angular_velocity >= self._torque_lut_omega[-1]: + return float(self._torque_lut[-1]) + return float( + np.interp(engine_angular_velocity, self._torque_lut_omega, self._torque_lut) + ) + + def _load_force(self, car_velocity: float) -> float: + return self._incline_force + self._drag_force_coeff * car_velocity * abs(car_velocity) - # Extract accelerations from breakdown - car_acceleration = system_breakdown.car.acceleration - engine_angular_accel = system_breakdown.engine.angular_acceleration - shift_acceleration = system_breakdown.cvt.acceleration + def _compute_dynamics( + self, y: list[float], full_shift: bool + ) -> tuple[float, float, float]: + car_velocity = y[0] + car_position = y[1] + shift_distance, shift_velocity = self._clamp_shift_state(y, full_shift) + engine_angular_velocity = y[4] + engine_angular_position = y[5] + + shift_distance_for_calc = min(shift_distance, self._pulley_calc_shift_max) + + cvt_ratio = self._lookup_cvt_ratio(shift_distance_for_calc) + cvt_ratio_derivative = self._lookup_cvt_ratio_rate( + shift_distance_for_calc, + shift_velocity, + ) + + engine_torque = self._lookup_engine_torque(engine_angular_velocity) + load_torque = self._load_force(car_velocity) * WHEEL_RADIUS + + wheel_angular_velocity = car_velocity / WHEEL_RADIUS + engine_to_wheel_ratio = cvt_ratio * GEARBOX_RATIO + engine_to_wheel_ratio_rate_of_change = cvt_ratio_derivative * GEARBOX_RATIO + + torque_demand = torque_demand_kernel( + engine_torque, + load_torque, + self._wheel_inertia, + wheel_angular_velocity, + engine_to_wheel_ratio, + engine_to_wheel_ratio_rate_of_change, + ENGINE_INERTIA, + ) + + state = self._scratch_state + state.car_velocity = car_velocity + state.car_position = car_position + state.shift_velocity = shift_velocity + state.shift_distance = shift_distance_for_calc + state.engine_angular_velocity = engine_angular_velocity + state.engine_angular_position = engine_angular_position + + t_max_prim = self._primary_pulley.calculate_max_torque(state) + t_max_sec = self._secondary_pulley.calculate_max_torque(state) + t_max_capacity = min(max(0.0, t_max_prim), max(0.0, t_max_sec)) + + relative_speed = slip_relative_speed_kernel( + engine_angular_velocity, + car_velocity * self._wheel_to_sec_ratio, + cvt_ratio, + ) + coupling_torque, _ = slip_coupling_torque_kernel( + relative_speed, + torque_demand, + t_max_capacity, + self._slip_speed_smoothing, + ) + + car_acceleration = ( + WHEEL_RADIUS * (coupling_torque * engine_to_wheel_ratio - load_torque) + ) / self._wheel_inertia + + engine_angular_accel = (engine_torque - coupling_torque) / ENGINE_INERTIA + + if full_shift: + return car_acceleration, engine_angular_accel, 0.0 + + primary_clamp, _ = self._primary_pulley.calculate_clamping_force(state) + _, _, primary_radial = self._primary_pulley.calculate_radial_force( + state, primary_clamp + ) + secondary_clamp, _ = self._secondary_pulley.calculate_clamping_force( + state, + torque=coupling_torque * cvt_ratio, + ) + _, _, secondary_radial = self._secondary_pulley.calculate_radial_force( + state, + secondary_clamp, + ) + + net_radial = primary_radial - secondary_radial + friction = self._cvt_shift_model._frictional_force(net_radial, shift_velocity) + shift_acceleration = (net_radial + friction) / self._cvt_moving_mass - # Prevent acceleration from pushing past boundaries (metal hitting metal) if shift_distance <= 0 and shift_acceleration < 0: - shift_acceleration = 0 + shift_acceleration = 0.0 elif shift_distance >= MAX_SHIFT and shift_acceleration > 0: - shift_acceleration = 0 + shift_acceleration = 0.0 + + return car_acceleration, engine_angular_accel, shift_acceleration + + def _evaluate_cvt_system(self, t: float, y: list[float]): + """Evaluate system dynamics (phase 1: not at full shift).""" + self._print_progress(t) + car_acceleration, engine_angular_accel, shift_acceleration = self._compute_dynamics( + y, + full_shift=False, + ) return [ car_acceleration, - state.car_velocity, + y[0], shift_acceleration, - state.shift_velocity, + y[2], engine_angular_accel, - state.engine_angular_velocity, + y[4], ] def _evaluate_full_shift_system(self, t: float, y: list[float]): """Evaluate system dynamics (phase 2: at full shift).""" - state = SystemState.from_array(y) self._print_progress(t) - # Force the shifting variables to remain constant at full shift. - state.shift_distance = MAX_SHIFT - state.shift_velocity = 0 - - # CRITICAL: Update the actual y array that scipy saves to CSV - constrained_y = state.to_array() - for i in range(len(y)): - y[i] = constrained_y[i] - - # Get system breakdown for full shift case - system_breakdown = self.system_model.get_breakdown(state) - - car_acceleration = system_breakdown.car.acceleration - engine_angular_accel = system_breakdown.engine.angular_acceleration + car_acceleration, engine_angular_accel, _ = self._compute_dynamics( + y, + full_shift=True, + ) return [ car_acceleration, - state.car_velocity, + y[0], 0, 0, engine_angular_accel, - state.engine_angular_velocity, + y[4], ] diff --git a/cvtModel/src/cvt_simulator/utils/numba_kernels.py b/cvtModel/src/cvt_simulator/utils/numba_kernels.py new file mode 100644 index 00000000..3c8e4000 --- /dev/null +++ b/cvtModel/src/cvt_simulator/utils/numba_kernels.py @@ -0,0 +1,171 @@ +import math + +import numpy as np + +try: + from numba import njit + + NUMBA_ENABLED = True + + def maybe_njit(*args, **kwargs): + return njit(*args, **kwargs) + +except ImportError: # pragma: no cover - exercised when numba is installed + NUMBA_ENABLED = False + + def maybe_njit(*args, **kwargs): + def decorator(func): + return func + + return decorator + + +@maybe_njit(cache=True, fastmath=True) +def slip_relative_speed_kernel( + primary_angular_velocity: float, + secondary_angular_velocity: float, + cvt_ratio: float, +) -> float: + return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) + + +@maybe_njit(cache=True, fastmath=True) +def slip_coupling_torque_kernel( + relative_speed: float, + torque_demand: float, + t_max_capacity: float, + slip_speed_smoothing: float, +) -> tuple[float, bool]: + coulomb_torque = t_max_capacity * math.tanh(relative_speed / slip_speed_smoothing) + alpha = min(max(abs(relative_speed) / slip_speed_smoothing, 0.0), 1.0) + torque_demand_clamped = min(max(torque_demand, -t_max_capacity), t_max_capacity) + coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque + return coupling_torque, alpha > 0.0 + + +@maybe_njit(cache=True, fastmath=True) +def torque_demand_kernel( + engine_torque: float, + load_torque: float, + wheel_inertia: float, + wheel_angular_velocity: float, + engine_to_wheel_ratio: float, + engine_to_wheel_ratio_rate_of_change: float, + engine_inertia: float, +) -> float: + eng_term = engine_torque * wheel_inertia + load_term = engine_inertia * load_torque * engine_to_wheel_ratio + shift_term = ( + engine_inertia + * wheel_inertia + * wheel_angular_velocity + * engine_to_wheel_ratio_rate_of_change + ) + + numerator = eng_term + load_term - shift_term + denominator = wheel_inertia + engine_inertia * engine_to_wheel_ratio**2 + return numerator / denominator + + +@maybe_njit(cache=True, fastmath=True) +def primary_flyweight_force_kernel( + flyweight_mass: float, + angular_velocity: float, + flyweight_radius: float, + ramp_slope: float, +) -> tuple[float, float, float, float]: + centrifugal_force = flyweight_mass * angular_velocity**2 * flyweight_radius + angle = math.atan(-ramp_slope) + angle_multiplier = math.tan(angle) + net = centrifugal_force * angle_multiplier + return angle, centrifugal_force, angle_multiplier, net + + +@maybe_njit(cache=True, fastmath=True) +def secondary_helix_force_kernel( + torque: float, + spring_torsion_torque: float, + secondary_radius: float, + ramp_slope: float, +) -> tuple[float, float, float]: + helix_angle = math.atan(-ramp_slope) + angle_multiplier = 2.0 * math.tan(helix_angle) * secondary_radius + net = (torque + spring_torsion_torque) / angle_multiplier + return helix_angle, angle_multiplier, net + + +@maybe_njit(cache=True, fastmath=True) +def radial_force_kernel( + clamping_force: float, + sheave_angle: float, + wrap_angle: float, + sec_angular_velocity: float, + sec_radius: float, + belt_cross_sectional_area: float, + rubber_density: float, +) -> tuple[float, float, float]: + radial_from_clamping = 2.0 * (clamping_force * math.tan(sheave_angle / 2.0)) / wrap_angle + + radial_from_centrifugal = ( + sec_angular_velocity**2 + * sec_radius**2 + * belt_cross_sectional_area + * rubber_density + ) + + total_radial = ( + 2.0 + * math.sin(wrap_angle / 2.0) + * (radial_from_clamping + radial_from_centrifugal) + ) + + return radial_from_clamping, radial_from_centrifugal, total_radial + + +@maybe_njit(cache=True, fastmath=True) +def max_torque_primary_kernel( + mu_effective: float, + wrap_angle: float, + total_radial: float, + radius: float, +) -> float: + exp_term = math.exp(mu_effective * wrap_angle) + capstan_term = (exp_term - 1.0) / (exp_term + 1.0) + radial_force_term = total_radial * radius / math.sin(wrap_angle / 2.0) + max_torque = capstan_term * radial_force_term + return max(0.0, max_torque) + + +@maybe_njit(cache=True, fastmath=True) +def max_torque_secondary_kernel( + spring_comp_force: float, + spring_tors_torque: float, + sheave_angle: float, + radial_from_centrifugal: float, + wrap_angle: float, + radius: float, + mu_effective: float, + cvt_ratio: float, + helix_radius: float, + helix_angle: float, +) -> float: + spring_force_term = (spring_comp_force + spring_tors_torque) * math.tan( + sheave_angle / 2.0 + ) + + centrifugal_force = radial_from_centrifugal * wrap_angle / 2.0 + + exp_term = math.exp(mu_effective * wrap_angle) + capstan_term = (wrap_angle / (4.0 * radius)) * (exp_term + 1.0) / (exp_term - 1.0) + + transmission_term = ( + 2.0 + * cvt_ratio + * (helix_radius * math.tan(helix_angle)) + * math.tan(sheave_angle / 2.0) + ) + + numerator = centrifugal_force + spring_force_term + denominator = capstan_term - transmission_term + max_torque = numerator / denominator + return max(0.0, max_torque) From 056871bee50fd7f76a394276c84d209cd83bc3f1 Mon Sep 17 00:00:00 2001 From: DO-Ui <68296986+DO-Ui@users.noreply.github.com> Date: Fri, 20 Feb 2026 13:52:07 -0500 Subject: [PATCH 2/3] In-place numba optimizations --- .../models/pulley/primary_pulley_flyweight.py | 37 +- .../models/pulley/pulley_interface.py | 43 ++- .../secondary_pulley_torque_reactive.py | 64 +++- .../src/cvt_simulator/models/slip_model.py | 72 +++- .../src/cvt_simulator/simulation_runner.py | 325 +++--------------- .../src/cvt_simulator/utils/numba_kernels.py | 171 --------- .../src/cvt_simulator/utils/numba_utils.py | 16 + 7 files changed, 256 insertions(+), 472 deletions(-) delete mode 100644 cvtModel/src/cvt_simulator/utils/numba_kernels.py create mode 100644 cvtModel/src/cvt_simulator/utils/numba_utils.py diff --git a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py index 03c90e08..b12f3705 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py +++ b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py @@ -18,10 +18,35 @@ INITIAL_FLYWEIGHT_RADIUS, ) from cvt_simulator.utils.system_state import SystemState -from cvt_simulator.utils.numba_kernels import ( - primary_flyweight_force_kernel, - max_torque_primary_kernel, -) +from cvt_simulator.utils.numba_utils import maybe_njit + + +@maybe_njit(cache=True, fastmath=True) +def _primary_flyweight_force_kernel( + flyweight_mass: float, + angular_velocity: float, + flyweight_radius: float, + ramp_slope: float, +) -> tuple[float, float, float, float]: + centrifugal_force = flyweight_mass * angular_velocity**2 * flyweight_radius + angle = math.atan(-ramp_slope) + angle_multiplier = math.tan(angle) + net = centrifugal_force * angle_multiplier + return angle, centrifugal_force, angle_multiplier, net + + +@maybe_njit(cache=True, fastmath=True) +def _max_torque_primary_kernel( + mu_effective: float, + wrap_angle: float, + total_radial: float, + radius: float, +) -> float: + exp_term = math.exp(mu_effective * wrap_angle) + capstan_term = (exp_term - 1.0) / (exp_term + 1.0) + radial_force_term = total_radial * radius / math.sin(wrap_angle / 2.0) + max_torque = capstan_term * radial_force_term + return max(0.0, max_torque) # TODO: Remove this code @@ -156,7 +181,7 @@ def calculate_max_torque( wrap_angle = self._get_wrap_angle(state.shift_distance) radius = self._get_radius(state.shift_distance) - return max_torque_primary_kernel(self.μ, wrap_angle, total_radial, radius) + return _max_torque_primary_kernel(self.μ, wrap_angle, total_radial, radius) # Private helper methods for force calculations @@ -174,7 +199,7 @@ def _calculate_flyweight_force( shift_distance ) - angle, centrifugal_force, angle_multiplier, net = primary_flyweight_force_kernel( + angle, centrifugal_force, angle_multiplier, net = _primary_flyweight_force_kernel( self.flyweight_mass, angular_velocity, flyweight_radius, diff --git a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py index 066ff584..0ea38b6c 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py +++ b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py @@ -35,7 +35,30 @@ ) from cvt_simulator.models.dataTypes import PulleyState, PulleyForces, PulleyBreakdowns from cvt_simulator.utils.theoretical_models import TheoreticalModels as tm -from cvt_simulator.utils.numba_kernels import radial_force_kernel +from cvt_simulator.utils.numba_utils import maybe_njit + + +@maybe_njit(cache=True, fastmath=True) +def _radial_force_kernel( + clamping_force: float, + sheave_angle: float, + wrap_angle: float, + sec_angular_velocity: float, + sec_radius: float, + belt_cross_sectional_area: float, + rubber_density: float, +) -> tuple[float, float, float]: + radial_from_clamping = 2.0 * (clamping_force * np.tan(sheave_angle / 2.0)) / wrap_angle + + radial_from_centrifugal = ( + sec_angular_velocity**2 * sec_radius**2 * belt_cross_sectional_area * rubber_density + ) + + total_radial = ( + 2.0 * np.sin(wrap_angle / 2.0) * (radial_from_clamping + radial_from_centrifugal) + ) + + return radial_from_clamping, radial_from_centrifugal, total_radial def get_kwarg(kwargs: dict[str, Any], key: str, default: Any = None) -> Any: @@ -183,16 +206,14 @@ def calculate_radial_force( wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS sec_angular_velocity = state.car_velocity * wheel_to_sec_ratio - radial_from_clamping, radial_from_centrifugal, total_radial = ( - radial_force_kernel( - clamping_force, - SHEAVE_ANGLE, - wrap_angle, - sec_angular_velocity, - sec_radius, - BELT_CROSS_SECTIONAL_AREA, - RUBBER_DENSITY, - ) + radial_from_clamping, radial_from_centrifugal, total_radial = _radial_force_kernel( + clamping_force, + SHEAVE_ANGLE, + wrap_angle, + sec_angular_velocity, + sec_radius, + BELT_CROSS_SECTIONAL_AREA, + RUBBER_DENSITY, ) return radial_from_clamping, radial_from_centrifugal, total_radial diff --git a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py index 333070c3..9b4dea23 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py +++ b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py @@ -16,10 +16,55 @@ ) from cvt_simulator.models.ramps import LinearSegment, PiecewiseRamp from cvt_simulator.utils.system_state import SystemState -from cvt_simulator.utils.numba_kernels import ( - secondary_helix_force_kernel, - max_torque_secondary_kernel, -) +from cvt_simulator.utils.numba_utils import maybe_njit + + +@maybe_njit(cache=True, fastmath=True) +def _secondary_helix_force_kernel( + torque: float, + spring_torsion_torque: float, + secondary_radius: float, + ramp_slope: float, +) -> tuple[float, float, float]: + helix_angle = np.arctan(-ramp_slope) + angle_multiplier = 2.0 * np.tan(helix_angle) * secondary_radius + net = (torque + spring_torsion_torque) / angle_multiplier + return helix_angle, angle_multiplier, net + + +@maybe_njit(cache=True, fastmath=True) +def _max_torque_secondary_kernel( + spring_comp_force: float, + spring_tors_torque: float, + sheave_angle: float, + radial_from_centrifugal: float, + wrap_angle: float, + radius: float, + mu_effective: float, + cvt_ratio: float, + helix_radius: float, + helix_angle: float, +) -> float: + spring_force_term = (spring_comp_force + spring_tors_torque) * np.tan( + sheave_angle / 2.0 + ) + + centrifugal_force = radial_from_centrifugal * wrap_angle / 2.0 + + exp_term = np.exp(mu_effective * wrap_angle) + capstan_term = (wrap_angle / (4.0 * radius)) * (exp_term + 1.0) / (exp_term - 1.0) + + transmission_term = ( + 2.0 + * cvt_ratio + * (helix_radius * np.tan(helix_angle)) + * np.tan(sheave_angle / 2.0) + ) + + numerator = centrifugal_force + spring_force_term + denominator = capstan_term - transmission_term + max_torque = numerator / denominator + return max(0.0, max_torque) # TODO: Remove this code @@ -87,6 +132,8 @@ def __init__( self.helix_radius = HELIX_RADIUS # This ramp needs to have a unique x for every height self.ramp = ramp + last_segment = self.ramp.segments[-1] + self._max_supported_shift = abs(self.ramp.height(last_segment.x_end)) def calculate_clamping_force( self, state: SystemState, **kwargs @@ -168,7 +215,7 @@ def calculate_max_torque( # Torque transmission term (feedback loop) cvt_ratio = tm.current_cvt_ratio(shift_distance) helix_angle = helix_breakdown.angle - return max_torque_secondary_kernel( + return _max_torque_secondary_kernel( spring_comp_force, spring_tors_torque, SHEAVE_ANGLE, @@ -193,7 +240,7 @@ def _calculate_helix_force( """ # Clamp shift distance to valid range # TODO: Remove - shift_distance = np.clip(shift_distance, 0, MAX_SHIFT) + shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) # Calculate torsion spring torque (resists cam rotation) spring_torque_breakdown = self._calculate_spring_tors_torque(shift_distance) @@ -201,7 +248,7 @@ def _calculate_helix_force( # Effective radius at current shift position secondary_radius = tm.outer_sec_radius(shift_distance) - BELT_HEIGHT / 2 - helix_angle, angle_multiplier, net = secondary_helix_force_kernel( + helix_angle, angle_multiplier, net = _secondary_helix_force_kernel( torque, spring_torque_breakdown.net, secondary_radius, @@ -238,7 +285,7 @@ def _calculate_spring_tors_torque( """Calculate torsion spring torque (resists helix cam rotation).""" # Clamp shift distance to valid range # TODO: Remove - shift_distance = np.clip(shift_distance, 0, MAX_SHIFT) + shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) # Calculate cam rotation from shift distance (approximation) # TODO: Improve relationship between shift and rotation @@ -266,6 +313,7 @@ def _calculate_rotation(self, shift_distance: float) -> float: Rotation angle [rad] """ # Find x position that corresponds to this height + shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) x_position = self.ramp.find_x_at_height(-shift_distance) # Get slope at that x position return x_position / HELIX_RADIUS diff --git a/cvtModel/src/cvt_simulator/models/slip_model.py b/cvtModel/src/cvt_simulator/models/slip_model.py index 7d72df82..421bf9ea 100644 --- a/cvtModel/src/cvt_simulator/models/slip_model.py +++ b/cvtModel/src/cvt_simulator/models/slip_model.py @@ -1,3 +1,4 @@ +import math import numpy as np from cvt_simulator.models.dataTypes import SlipBreakdown from cvt_simulator.models.external_load_model import LoadModel @@ -16,11 +17,54 @@ from cvt_simulator.constants.constants import ( RUBBER_ALUMINUM_STATIC_FRICTION, ) -from cvt_simulator.utils.numba_kernels import ( - slip_relative_speed_kernel, - slip_coupling_torque_kernel, - torque_demand_kernel, -) +from cvt_simulator.utils.numba_utils import maybe_njit + + +@maybe_njit(cache=True, fastmath=True) +def _relative_speed_kernel( + primary_angular_velocity: float, + secondary_angular_velocity: float, + cvt_ratio: float, +) -> float: + return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) + + +@maybe_njit(cache=True, fastmath=True) +def _coupling_torque_kernel( + relative_speed: float, + torque_demand: float, + t_max_capacity: float, + slip_speed_smoothing: float, +) -> tuple[float, bool]: + coulomb_torque = t_max_capacity * math.tanh(relative_speed / slip_speed_smoothing) + alpha = min(max(abs(relative_speed) / slip_speed_smoothing, 0.0), 1.0) + torque_demand_clamped = min(max(torque_demand, -t_max_capacity), t_max_capacity) + coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque + return coupling_torque, alpha > 0.0 + + +@maybe_njit(cache=True, fastmath=True) +def _torque_demand_kernel( + engine_torque: float, + load_torque: float, + wheel_inertia: float, + wheel_angular_velocity: float, + engine_to_wheel_ratio: float, + engine_to_wheel_ratio_rate_of_change: float, + engine_inertia: float, +) -> float: + eng_term = engine_torque * wheel_inertia + load_term = engine_inertia * load_torque * engine_to_wheel_ratio + shift_term = ( + engine_inertia + * wheel_inertia + * wheel_angular_velocity + * engine_to_wheel_ratio_rate_of_change + ) + + numerator = eng_term + load_term - shift_term + denominator = wheel_inertia + engine_inertia * engine_to_wheel_ratio**2 + return numerator / denominator class SlipModel: @@ -62,12 +106,12 @@ def get_breakdown(self, state: SystemState) -> SlipBreakdown: torque_demand = self.get_torque_demand(state) wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS - relative_speed = slip_relative_speed_kernel( + relative_speed = self._relative_speed( state.engine_angular_velocity, state.car_velocity * wheel_to_sec_ratio, tm.current_cvt_ratio(state.shift_distance), ) - coupling_torque, _ = slip_coupling_torque_kernel( + coupling_torque, _ = _coupling_torque_kernel( relative_speed, torque_demand, t_max_capacity, @@ -105,7 +149,7 @@ def get_torque_demand(self, state: SystemState): * GEARBOX_RATIO ) - coupling_torque = torque_demand_kernel( + coupling_torque = _torque_demand_kernel( engine_torque, load_torque, wheel_inertia, @@ -141,6 +185,18 @@ def calculate_t_max(self, state: SystemState) -> tuple[float, float]: secondary_t_max = max(0, secondary_t_max) return primary_t_max, secondary_t_max + def _relative_speed( + self, + primary_angular_velocity: float, + secondary_angular_velocity: float, + cvt_ratio: float, + ) -> float: + return _relative_speed_kernel( + primary_angular_velocity, + secondary_angular_velocity, + cvt_ratio, + ) + # def _is_slipping( # self, # primary_angular_velocity: float, diff --git a/cvtModel/src/cvt_simulator/simulation_runner.py b/cvtModel/src/cvt_simulator/simulation_runner.py index ac202cf0..ea284c42 100644 --- a/cvtModel/src/cvt_simulator/simulation_runner.py +++ b/cvtModel/src/cvt_simulator/simulation_runner.py @@ -1,5 +1,4 @@ import sys -import math import numpy as np from typing import Callable, Optional from scipy.integrate import solve_ivp @@ -10,16 +9,14 @@ GEARBOX_RATIO, WHEEL_RADIUS, MAX_SHIFT, - ENGINE_INERTIA, - DRIVELINE_INERTIA, ) from cvt_simulator.utils.conversions import rpm_to_rad_s from cvt_simulator.utils.theoretical_models import TheoreticalModels as tm -from cvt_simulator.constants.engine_specs import safe_torque_curve -from cvt_simulator.utils.numba_kernels import ( - slip_relative_speed_kernel, - slip_coupling_torque_kernel, - torque_demand_kernel, +from cvt_simulator.utils.simulation_constraints import ( + car_velocity_constraint_event, + get_shift_steady_event, + get_back_shift_event, + shift_constraint_event, ) @@ -55,63 +52,6 @@ def __init__( self.progress_callback = progress_callback self._last_callback_percent = -1.0 - # Cached model references for faster ODE RHS evaluation - self._slip_model = system_model.slip_model - self._cvt_shift_model = system_model.cvt_shift_model - self._primary_pulley = self._slip_model.primary_pulley - self._secondary_pulley = self._slip_model.secondary_pulley - self._engine_model = self._slip_model.engine_model - self._load_model = self._slip_model.load_model - - self._wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS - self._wheel_inertia = DRIVELINE_INERTIA + self._slip_model.car_mass * ( - WHEEL_RADIUS**2 - ) - self._cvt_moving_mass = self._cvt_shift_model.cvt_moving_mass - self._slip_speed_smoothing = self._slip_model.slip_speed_smoothing - - self._incline_force = ( - self._load_model.car_mass - * self._load_model.g - * math.sin(self._load_model.incline_angle) - ) - self._drag_force_coeff = ( - 0.5 - * self._load_model.air_density - * self._load_model.frontal_area - * self._load_model.drag_coefficient - ) - - # Precompute CVT ratio and derivative lookup tables to avoid per-step root finding - self._ratio_lut_shift = np.linspace(0.0, MAX_SHIFT, 1024) - self._ratio_lut = np.array( - [tm.current_cvt_ratio(d) for d in self._ratio_lut_shift], dtype=float - ) - self._ratio_rate_per_vel_lut = np.array( - [ - tm.current_cvt_ratio_rate_of_change(d, 1.0) - for d in self._ratio_lut_shift - ], - dtype=float, - ) - - # Precompute torque lookup table to avoid repeated scipy interp1d overhead - self._torque_lut_omega = np.linspace(0.0, rpm_to_rad_s(6000), 2048) - self._torque_lut = np.asarray(safe_torque_curve(self._torque_lut_omega), dtype=float) - - # Reused mutable state object to minimize allocation churn - self._scratch_state = SystemState() - - # Secondary ramp can have a slightly smaller usable max than MAX_SHIFT. - # Clamp internal pulley calculations to avoid inverse-height edge errors. - secondary_last_segment = self._secondary_pulley.ramp.segments[-1] - secondary_max_shift = abs( - self._secondary_pulley.ramp.height(secondary_last_segment.x_end) - ) - self._pulley_calc_shift_max = min(MAX_SHIFT, secondary_max_shift) - - - def run_simulation(self) -> SimulationResult: """Run the simulation and return results.""" cvt_system_ode = self._get_ode_function() @@ -127,9 +67,9 @@ def run_simulation(self) -> SimulationResult: # Phase 1: Normal shifting until full shift is reached events_phase1 = [ - self._get_shift_steady_event(), - self._get_car_velocity_constraint_event(), - self._get_shift_constraint_event(), + get_shift_steady_event(self.system_model), + car_velocity_constraint_event, + shift_constraint_event, ] time_eval_phase1 = time_eval[time_eval >= current_time] @@ -162,8 +102,8 @@ def run_simulation(self) -> SimulationResult: # At full shift, check for back-shift event events_phase2 = [ - self._get_back_shift_event(), - self._get_car_velocity_constraint_event(), + get_back_shift_event(self.system_model), + car_velocity_constraint_event, ] solution_phase2 = self._solve( @@ -189,9 +129,9 @@ def run_simulation(self) -> SimulationResult: break events_phase3 = [ - self._get_shift_steady_event(), - self._get_car_velocity_constraint_event(), - self._get_shift_constraint_event(), + get_shift_steady_event(self.system_model), + car_velocity_constraint_event, + shift_constraint_event, ] solution_phase3 = self._solve( @@ -273,224 +213,73 @@ def _print_progress(self, t): self._last_callback_percent = rounded_percent self.progress_callback(progress_percent) - def _get_shift_constraint_event(self): - def shift_constraint_event(t, y): - shift_velocity = y[2] - shift_distance = y[3] - - if shift_distance < 0: - y[3] = 0.0 - y[2] = max(0.0, shift_velocity) - elif shift_distance > MAX_SHIFT: - y[3] = MAX_SHIFT - y[2] = min(0.0, shift_velocity) - - return 1.0 - - return shift_constraint_event - - def _get_car_velocity_constraint_event(self): - def car_velocity_event(t, y): - return y[0] - - car_velocity_event.terminal = True - car_velocity_event.direction = -1 - return car_velocity_event - - def _get_shift_steady_event(self): - def shift_steady_event(t, y): - tol = 1e-5 - if y[3] < MAX_SHIFT - tol: - return -tol - - shift_velocity = y[2] - shift_distance = y[3] - if shift_distance < 0: - y[3] = 0.0 - y[2] = max(0.0, shift_velocity) - elif shift_distance > MAX_SHIFT: - y[3] = MAX_SHIFT - y[2] = min(0.0, shift_velocity) - - _, _, shift_acceleration = self._compute_dynamics(y, full_shift=False) - return shift_acceleration - - shift_steady_event.terminal = True - shift_steady_event.direction = 1 - return shift_steady_event - - def _get_back_shift_event(self): - def back_shift_event(t, y): - if y[3] < MAX_SHIFT - 1e-5: - return 1.0 - - _, _, shift_acceleration = self._compute_dynamics(y, full_shift=False) - return shift_acceleration + 5.0 - - back_shift_event.terminal = True - back_shift_event.direction = -1 - return back_shift_event - - def _clamp_shift_state(self, y: list[float], full_shift: bool) -> tuple[float, float]: - if full_shift: - y[3] = MAX_SHIFT - y[2] = 0.0 - return MAX_SHIFT, 0.0 - - shift_velocity = y[2] - shift_distance = y[3] + def _evaluate_cvt_system(self, t: float, y: list[float]): + """Evaluate system dynamics (phase 1: not at full shift).""" + state = SystemState.from_array(y) + self._print_progress(t) + # TODO: Remove this (should be handled by constraints) + shift_velocity = state.shift_velocity + shift_distance = state.shift_distance if shift_distance <= 0: - y[3] = 0.0 - y[2] = max(0.0, shift_velocity) - elif shift_distance > MAX_SHIFT: - y[3] = MAX_SHIFT - y[2] = min(0.0, shift_velocity) - - return y[3], y[2] - - def _lookup_cvt_ratio(self, shift_distance: float) -> float: - return float(np.interp(shift_distance, self._ratio_lut_shift, self._ratio_lut)) - - def _lookup_cvt_ratio_rate(self, shift_distance: float, shift_velocity: float) -> float: - rate_per_velocity = float( - np.interp(shift_distance, self._ratio_lut_shift, self._ratio_rate_per_vel_lut) - ) - return rate_per_velocity * shift_velocity - - def _lookup_engine_torque(self, engine_angular_velocity: float) -> float: - if engine_angular_velocity <= self._torque_lut_omega[0]: - return float(self._torque_lut[0]) - if engine_angular_velocity >= self._torque_lut_omega[-1]: - return float(self._torque_lut[-1]) - return float( - np.interp(engine_angular_velocity, self._torque_lut_omega, self._torque_lut) - ) - - def _load_force(self, car_velocity: float) -> float: - return self._incline_force + self._drag_force_coeff * car_velocity * abs(car_velocity) - - def _compute_dynamics( - self, y: list[float], full_shift: bool - ) -> tuple[float, float, float]: - car_velocity = y[0] - car_position = y[1] - shift_distance, shift_velocity = self._clamp_shift_state(y, full_shift) - engine_angular_velocity = y[4] - engine_angular_position = y[5] + state.shift_distance = 0 + state.shift_velocity = max(0, shift_velocity) - shift_distance_for_calc = min(shift_distance, self._pulley_calc_shift_max) - - cvt_ratio = self._lookup_cvt_ratio(shift_distance_for_calc) - cvt_ratio_derivative = self._lookup_cvt_ratio_rate( - shift_distance_for_calc, - shift_velocity, - ) - - engine_torque = self._lookup_engine_torque(engine_angular_velocity) - load_torque = self._load_force(car_velocity) * WHEEL_RADIUS - - wheel_angular_velocity = car_velocity / WHEEL_RADIUS - engine_to_wheel_ratio = cvt_ratio * GEARBOX_RATIO - engine_to_wheel_ratio_rate_of_change = cvt_ratio_derivative * GEARBOX_RATIO - - torque_demand = torque_demand_kernel( - engine_torque, - load_torque, - self._wheel_inertia, - wheel_angular_velocity, - engine_to_wheel_ratio, - engine_to_wheel_ratio_rate_of_change, - ENGINE_INERTIA, - ) - - state = self._scratch_state - state.car_velocity = car_velocity - state.car_position = car_position - state.shift_velocity = shift_velocity - state.shift_distance = shift_distance_for_calc - state.engine_angular_velocity = engine_angular_velocity - state.engine_angular_position = engine_angular_position - - t_max_prim = self._primary_pulley.calculate_max_torque(state) - t_max_sec = self._secondary_pulley.calculate_max_torque(state) - t_max_capacity = min(max(0.0, t_max_prim), max(0.0, t_max_sec)) - - relative_speed = slip_relative_speed_kernel( - engine_angular_velocity, - car_velocity * self._wheel_to_sec_ratio, - cvt_ratio, - ) - coupling_torque, _ = slip_coupling_torque_kernel( - relative_speed, - torque_demand, - t_max_capacity, - self._slip_speed_smoothing, - ) - - car_acceleration = ( - WHEEL_RADIUS * (coupling_torque * engine_to_wheel_ratio - load_torque) - ) / self._wheel_inertia + elif shift_distance > MAX_SHIFT: + state.shift_distance = MAX_SHIFT + state.shift_velocity = min(0, shift_velocity) - engine_angular_accel = (engine_torque - coupling_torque) / ENGINE_INERTIA + constrained_y = state.to_array() + for i in range(len(y)): + y[i] = constrained_y[i] - if full_shift: - return car_acceleration, engine_angular_accel, 0.0 + # Get system breakdown (this calculates everything in correct order) + system_breakdown = self.system_model.get_breakdown(state) - primary_clamp, _ = self._primary_pulley.calculate_clamping_force(state) - _, _, primary_radial = self._primary_pulley.calculate_radial_force( - state, primary_clamp - ) - secondary_clamp, _ = self._secondary_pulley.calculate_clamping_force( - state, - torque=coupling_torque * cvt_ratio, - ) - _, _, secondary_radial = self._secondary_pulley.calculate_radial_force( - state, - secondary_clamp, - ) - - net_radial = primary_radial - secondary_radial - friction = self._cvt_shift_model._frictional_force(net_radial, shift_velocity) - shift_acceleration = (net_radial + friction) / self._cvt_moving_mass + # Extract accelerations from breakdown + car_acceleration = system_breakdown.car.acceleration + engine_angular_accel = system_breakdown.engine.angular_acceleration + shift_acceleration = system_breakdown.cvt.acceleration + # Prevent acceleration from pushing past boundaries (metal hitting metal) if shift_distance <= 0 and shift_acceleration < 0: - shift_acceleration = 0.0 + shift_acceleration = 0 elif shift_distance >= MAX_SHIFT and shift_acceleration > 0: - shift_acceleration = 0.0 - - return car_acceleration, engine_angular_accel, shift_acceleration - - def _evaluate_cvt_system(self, t: float, y: list[float]): - """Evaluate system dynamics (phase 1: not at full shift).""" - self._print_progress(t) - car_acceleration, engine_angular_accel, shift_acceleration = self._compute_dynamics( - y, - full_shift=False, - ) + shift_acceleration = 0 return [ car_acceleration, - y[0], + state.car_velocity, shift_acceleration, - y[2], + state.shift_velocity, engine_angular_accel, - y[4], + state.engine_angular_velocity, ] def _evaluate_full_shift_system(self, t: float, y: list[float]): """Evaluate system dynamics (phase 2: at full shift).""" + state = SystemState.from_array(y) self._print_progress(t) - car_acceleration, engine_angular_accel, _ = self._compute_dynamics( - y, - full_shift=True, - ) + # Force the shifting variables to remain constant at full shift. + state.shift_distance = MAX_SHIFT + state.shift_velocity = 0 + + # CRITICAL: Update the actual y array that scipy saves to CSV + constrained_y = state.to_array() + for i in range(len(y)): + y[i] = constrained_y[i] + + # Get system breakdown for full shift case + system_breakdown = self.system_model.get_breakdown(state) + + car_acceleration = system_breakdown.car.acceleration + engine_angular_accel = system_breakdown.engine.angular_acceleration return [ car_acceleration, - y[0], + state.car_velocity, 0, 0, engine_angular_accel, - y[4], + state.engine_angular_velocity, ] diff --git a/cvtModel/src/cvt_simulator/utils/numba_kernels.py b/cvtModel/src/cvt_simulator/utils/numba_kernels.py deleted file mode 100644 index 3c8e4000..00000000 --- a/cvtModel/src/cvt_simulator/utils/numba_kernels.py +++ /dev/null @@ -1,171 +0,0 @@ -import math - -import numpy as np - -try: - from numba import njit - - NUMBA_ENABLED = True - - def maybe_njit(*args, **kwargs): - return njit(*args, **kwargs) - -except ImportError: # pragma: no cover - exercised when numba is installed - NUMBA_ENABLED = False - - def maybe_njit(*args, **kwargs): - def decorator(func): - return func - - return decorator - - -@maybe_njit(cache=True, fastmath=True) -def slip_relative_speed_kernel( - primary_angular_velocity: float, - secondary_angular_velocity: float, - cvt_ratio: float, -) -> float: - return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) - - -@maybe_njit(cache=True, fastmath=True) -def slip_coupling_torque_kernel( - relative_speed: float, - torque_demand: float, - t_max_capacity: float, - slip_speed_smoothing: float, -) -> tuple[float, bool]: - coulomb_torque = t_max_capacity * math.tanh(relative_speed / slip_speed_smoothing) - alpha = min(max(abs(relative_speed) / slip_speed_smoothing, 0.0), 1.0) - torque_demand_clamped = min(max(torque_demand, -t_max_capacity), t_max_capacity) - coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque - return coupling_torque, alpha > 0.0 - - -@maybe_njit(cache=True, fastmath=True) -def torque_demand_kernel( - engine_torque: float, - load_torque: float, - wheel_inertia: float, - wheel_angular_velocity: float, - engine_to_wheel_ratio: float, - engine_to_wheel_ratio_rate_of_change: float, - engine_inertia: float, -) -> float: - eng_term = engine_torque * wheel_inertia - load_term = engine_inertia * load_torque * engine_to_wheel_ratio - shift_term = ( - engine_inertia - * wheel_inertia - * wheel_angular_velocity - * engine_to_wheel_ratio_rate_of_change - ) - - numerator = eng_term + load_term - shift_term - denominator = wheel_inertia + engine_inertia * engine_to_wheel_ratio**2 - return numerator / denominator - - -@maybe_njit(cache=True, fastmath=True) -def primary_flyweight_force_kernel( - flyweight_mass: float, - angular_velocity: float, - flyweight_radius: float, - ramp_slope: float, -) -> tuple[float, float, float, float]: - centrifugal_force = flyweight_mass * angular_velocity**2 * flyweight_radius - angle = math.atan(-ramp_slope) - angle_multiplier = math.tan(angle) - net = centrifugal_force * angle_multiplier - return angle, centrifugal_force, angle_multiplier, net - - -@maybe_njit(cache=True, fastmath=True) -def secondary_helix_force_kernel( - torque: float, - spring_torsion_torque: float, - secondary_radius: float, - ramp_slope: float, -) -> tuple[float, float, float]: - helix_angle = math.atan(-ramp_slope) - angle_multiplier = 2.0 * math.tan(helix_angle) * secondary_radius - net = (torque + spring_torsion_torque) / angle_multiplier - return helix_angle, angle_multiplier, net - - -@maybe_njit(cache=True, fastmath=True) -def radial_force_kernel( - clamping_force: float, - sheave_angle: float, - wrap_angle: float, - sec_angular_velocity: float, - sec_radius: float, - belt_cross_sectional_area: float, - rubber_density: float, -) -> tuple[float, float, float]: - radial_from_clamping = 2.0 * (clamping_force * math.tan(sheave_angle / 2.0)) / wrap_angle - - radial_from_centrifugal = ( - sec_angular_velocity**2 - * sec_radius**2 - * belt_cross_sectional_area - * rubber_density - ) - - total_radial = ( - 2.0 - * math.sin(wrap_angle / 2.0) - * (radial_from_clamping + radial_from_centrifugal) - ) - - return radial_from_clamping, radial_from_centrifugal, total_radial - - -@maybe_njit(cache=True, fastmath=True) -def max_torque_primary_kernel( - mu_effective: float, - wrap_angle: float, - total_radial: float, - radius: float, -) -> float: - exp_term = math.exp(mu_effective * wrap_angle) - capstan_term = (exp_term - 1.0) / (exp_term + 1.0) - radial_force_term = total_radial * radius / math.sin(wrap_angle / 2.0) - max_torque = capstan_term * radial_force_term - return max(0.0, max_torque) - - -@maybe_njit(cache=True, fastmath=True) -def max_torque_secondary_kernel( - spring_comp_force: float, - spring_tors_torque: float, - sheave_angle: float, - radial_from_centrifugal: float, - wrap_angle: float, - radius: float, - mu_effective: float, - cvt_ratio: float, - helix_radius: float, - helix_angle: float, -) -> float: - spring_force_term = (spring_comp_force + spring_tors_torque) * math.tan( - sheave_angle / 2.0 - ) - - centrifugal_force = radial_from_centrifugal * wrap_angle / 2.0 - - exp_term = math.exp(mu_effective * wrap_angle) - capstan_term = (wrap_angle / (4.0 * radius)) * (exp_term + 1.0) / (exp_term - 1.0) - - transmission_term = ( - 2.0 - * cvt_ratio - * (helix_radius * math.tan(helix_angle)) - * math.tan(sheave_angle / 2.0) - ) - - numerator = centrifugal_force + spring_force_term - denominator = capstan_term - transmission_term - max_torque = numerator / denominator - return max(0.0, max_torque) diff --git a/cvtModel/src/cvt_simulator/utils/numba_utils.py b/cvtModel/src/cvt_simulator/utils/numba_utils.py new file mode 100644 index 00000000..6da40a38 --- /dev/null +++ b/cvtModel/src/cvt_simulator/utils/numba_utils.py @@ -0,0 +1,16 @@ +try: + from numba import njit + + NUMBA_ENABLED = True + + def maybe_njit(*args, **kwargs): + return njit(*args, **kwargs) + +except ImportError: # pragma: no cover - exercised when numba is installed + NUMBA_ENABLED = False + + def maybe_njit(*args, **kwargs): + def decorator(func): + return func + + return decorator From 949d6332c2421ef502351b6832ebd9fed464d3b9 Mon Sep 17 00:00:00 2001 From: DO-Ui <68296986+DO-Ui@users.noreply.github.com> Date: Fri, 20 Feb 2026 18:03:15 -0500 Subject: [PATCH 3/3] Annotations Only Version --- cvtModel/README.md | 4 - cvtModel/setup.py | 6 +- .../models/pulley/primary_pulley_flyweight.py | 53 ++++----- .../models/pulley/pulley_interface.py | 49 +++----- .../secondary_pulley_torque_reactive.py | 105 ++++++------------ .../src/cvt_simulator/models/slip_model.py | 99 ++++++----------- 6 files changed, 105 insertions(+), 211 deletions(-) diff --git a/cvtModel/README.md b/cvtModel/README.md index 9eff9ccc..bcc26ce9 100644 --- a/cvtModel/README.md +++ b/cvtModel/README.md @@ -9,10 +9,6 @@ venv\Scripts\activate # or 'source venv/bin/activate' on Mac ```bash pip install -e .[dev] ``` -For optional runtime acceleration kernels: -```bash -pip install -e .[dev,numba] -``` 3. Run the tests ```bash coverage run -m unittest discover -s test/simulations -s test/utils diff --git a/cvtModel/setup.py b/cvtModel/setup.py index 83c9baeb..73d06f05 100644 --- a/cvtModel/setup.py +++ b/cvtModel/setup.py @@ -10,7 +10,8 @@ "scipy", "matplotlib", "pandas", - "pydantic" + "pydantic", + "numba" ], extras_require={ "dev": [ @@ -19,8 +20,5 @@ "pytest", "coverage" ], - "numba": [ - "numba" - ] }, ) diff --git a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py index b12f3705..15e3e733 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py +++ b/cvtModel/src/cvt_simulator/models/pulley/primary_pulley_flyweight.py @@ -18,35 +18,6 @@ INITIAL_FLYWEIGHT_RADIUS, ) from cvt_simulator.utils.system_state import SystemState -from cvt_simulator.utils.numba_utils import maybe_njit - - -@maybe_njit(cache=True, fastmath=True) -def _primary_flyweight_force_kernel( - flyweight_mass: float, - angular_velocity: float, - flyweight_radius: float, - ramp_slope: float, -) -> tuple[float, float, float, float]: - centrifugal_force = flyweight_mass * angular_velocity**2 * flyweight_radius - angle = math.atan(-ramp_slope) - angle_multiplier = math.tan(angle) - net = centrifugal_force * angle_multiplier - return angle, centrifugal_force, angle_multiplier, net - - -@maybe_njit(cache=True, fastmath=True) -def _max_torque_primary_kernel( - mu_effective: float, - wrap_angle: float, - total_radial: float, - radius: float, -) -> float: - exp_term = math.exp(mu_effective * wrap_angle) - capstan_term = (exp_term - 1.0) / (exp_term + 1.0) - radial_force_term = total_radial * radius / math.sin(wrap_angle / 2.0) - max_torque = capstan_term * radial_force_term - return max(0.0, max_torque) # TODO: Remove this code @@ -181,7 +152,14 @@ def calculate_max_torque( wrap_angle = self._get_wrap_angle(state.shift_distance) radius = self._get_radius(state.shift_distance) - return _max_torque_primary_kernel(self.μ, wrap_angle, total_radial, radius) + # Capstan equation with V-belt friction enhancement + exp_term = math.exp(self.μ * wrap_angle) + capstan_term = (exp_term - 1) / (exp_term + 1) + radial_force_term = total_radial * radius / np.sin(wrap_angle / 2) + + max_torque = capstan_term * radial_force_term + + return max(0.0, max_torque) # Ensure non-negative # Private helper methods for force calculations @@ -199,19 +177,28 @@ def _calculate_flyweight_force( shift_distance ) - angle, centrifugal_force, angle_multiplier, net = _primary_flyweight_force_kernel( + # Centrifugal force on flyweight: F = m * ω² * r + centrifugal_force = tm.centrifugal_force( self.flyweight_mass, angular_velocity, flyweight_radius, - self.ramp.slope(shift_distance), ) + # Ramp angle at current position + # Ramp is default negative slope, so negate for angle + # If a positive slope ramp is passed, it will generate a force against shifting + # which is expected for such a stupid ramp design. + angle = np.arctan(-self.ramp.slope(shift_distance)) + + # Convert centrifugal force to axial clamping force through ramp angle + net = centrifugal_force * np.tan(angle) + return flyweightForceBreakdown( radius=flyweight_radius, angular_velocity=angular_velocity, angle=angle, centrifugal_force=centrifugal_force, - angle_multiplier=angle_multiplier, + angle_multiplier=np.tan(angle), net=net, ) diff --git a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py index 0ea38b6c..6b0b0ca5 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py +++ b/cvtModel/src/cvt_simulator/models/pulley/pulley_interface.py @@ -35,30 +35,6 @@ ) from cvt_simulator.models.dataTypes import PulleyState, PulleyForces, PulleyBreakdowns from cvt_simulator.utils.theoretical_models import TheoreticalModels as tm -from cvt_simulator.utils.numba_utils import maybe_njit - - -@maybe_njit(cache=True, fastmath=True) -def _radial_force_kernel( - clamping_force: float, - sheave_angle: float, - wrap_angle: float, - sec_angular_velocity: float, - sec_radius: float, - belt_cross_sectional_area: float, - rubber_density: float, -) -> tuple[float, float, float]: - radial_from_clamping = 2.0 * (clamping_force * np.tan(sheave_angle / 2.0)) / wrap_angle - - radial_from_centrifugal = ( - sec_angular_velocity**2 * sec_radius**2 * belt_cross_sectional_area * rubber_density - ) - - total_radial = ( - 2.0 * np.sin(wrap_angle / 2.0) * (radial_from_clamping + radial_from_centrifugal) - ) - - return radial_from_clamping, radial_from_centrifugal, total_radial def get_kwarg(kwargs: dict[str, Any], key: str, default: Any = None) -> Any: @@ -196,6 +172,11 @@ def calculate_radial_force( # radius = self._get_radius(state.shift_distance) # angular_velocity = self._get_angular_velocity(state) + # Radial force from pulley clamping (through V-belt wedging) + radial_from_clamping = ( + 2 * (clamping_force * np.tan(SHEAVE_ANGLE / 2)) / wrap_angle + ) + # To future Kai: # Trust me! If you're using your current assumption # (Which the belt isn't slipping from the secondary, i.e. belt speed = sec speed) @@ -206,14 +187,18 @@ def calculate_radial_force( wheel_to_sec_ratio = GEARBOX_RATIO / WHEEL_RADIUS sec_angular_velocity = state.car_velocity * wheel_to_sec_ratio - radial_from_clamping, radial_from_centrifugal, total_radial = _radial_force_kernel( - clamping_force, - SHEAVE_ANGLE, - wrap_angle, - sec_angular_velocity, - sec_radius, - BELT_CROSS_SECTIONAL_AREA, - RUBBER_DENSITY, + radial_from_centrifugal = ( + sec_angular_velocity**2 + * sec_radius**2 + * BELT_CROSS_SECTIONAL_AREA + * RUBBER_DENSITY + ) + + # Total radial force (determines friction capacity) + total_radial = ( + 2 + * np.sin(wrap_angle / 2) + * (radial_from_clamping + radial_from_centrifugal) ) return radial_from_clamping, radial_from_centrifugal, total_radial diff --git a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py index 9b4dea23..9c4faab8 100644 --- a/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py +++ b/cvtModel/src/cvt_simulator/models/pulley/secondary_pulley_torque_reactive.py @@ -16,55 +16,6 @@ ) from cvt_simulator.models.ramps import LinearSegment, PiecewiseRamp from cvt_simulator.utils.system_state import SystemState -from cvt_simulator.utils.numba_utils import maybe_njit - - -@maybe_njit(cache=True, fastmath=True) -def _secondary_helix_force_kernel( - torque: float, - spring_torsion_torque: float, - secondary_radius: float, - ramp_slope: float, -) -> tuple[float, float, float]: - helix_angle = np.arctan(-ramp_slope) - angle_multiplier = 2.0 * np.tan(helix_angle) * secondary_radius - net = (torque + spring_torsion_torque) / angle_multiplier - return helix_angle, angle_multiplier, net - - -@maybe_njit(cache=True, fastmath=True) -def _max_torque_secondary_kernel( - spring_comp_force: float, - spring_tors_torque: float, - sheave_angle: float, - radial_from_centrifugal: float, - wrap_angle: float, - radius: float, - mu_effective: float, - cvt_ratio: float, - helix_radius: float, - helix_angle: float, -) -> float: - spring_force_term = (spring_comp_force + spring_tors_torque) * np.tan( - sheave_angle / 2.0 - ) - - centrifugal_force = radial_from_centrifugal * wrap_angle / 2.0 - - exp_term = np.exp(mu_effective * wrap_angle) - capstan_term = (wrap_angle / (4.0 * radius)) * (exp_term + 1.0) / (exp_term - 1.0) - - transmission_term = ( - 2.0 - * cvt_ratio - * (helix_radius * np.tan(helix_angle)) - * np.tan(sheave_angle / 2.0) - ) - - numerator = centrifugal_force + spring_force_term - denominator = capstan_term - transmission_term - max_torque = numerator / denominator - return max(0.0, max_torque) # TODO: Remove this code @@ -132,8 +83,6 @@ def __init__( self.helix_radius = HELIX_RADIUS # This ramp needs to have a unique x for every height self.ramp = ramp - last_segment = self.ramp.segments[-1] - self._max_supported_shift = abs(self.ramp.height(last_segment.x_end)) def calculate_clamping_force( self, state: SystemState, **kwargs @@ -209,25 +158,36 @@ def calculate_max_torque( helix_breakdown = self._calculate_helix_force(0, shift_distance) spring_tors_torque = helix_breakdown.net + # Convert to radial force contribution + spring_force_term = (spring_comp_force + spring_tors_torque) * np.tan( + SHEAVE_ANGLE / 2 + ) + # Centrifugal force contribution, used built in calc with 0 clamp (since we only need centrifugal) _, radial_from_centrifugal, _ = self.calculate_radial_force(state, 0) + centrifugal_force = radial_from_centrifugal * wrap_angle / 2 + + # Capstan term + exp_term = np.exp(self.μ * wrap_angle) + capstan_term = (wrap_angle / (4 * radius)) * (exp_term + 1) / (exp_term - 1) # Torque transmission term (feedback loop) cvt_ratio = tm.current_cvt_ratio(shift_distance) helix_angle = helix_breakdown.angle - return _max_torque_secondary_kernel( - spring_comp_force, - spring_tors_torque, - SHEAVE_ANGLE, - radial_from_centrifugal, - wrap_angle, - radius, - self.μ, - cvt_ratio, - HELIX_RADIUS, - helix_angle, + transmission_term = ( + 2 + * cvt_ratio + * (HELIX_RADIUS * np.tan(helix_angle)) + * np.tan(SHEAVE_ANGLE / 2) ) + # Solve for max torque (equilibrium of torque feedback loop) + numerator = centrifugal_force + spring_force_term + denominator = capstan_term - transmission_term + max_torque = numerator / denominator + + return max(0.0, max_torque) # Ensure non-negative + # Private helper methods for force calculations def _calculate_helix_force( @@ -240,7 +200,7 @@ def _calculate_helix_force( """ # Clamp shift distance to valid range # TODO: Remove - shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) + shift_distance = np.clip(shift_distance, 0, MAX_SHIFT) # Calculate torsion spring torque (resists cam rotation) spring_torque_breakdown = self._calculate_spring_tors_torque(shift_distance) @@ -248,12 +208,16 @@ def _calculate_helix_force( # Effective radius at current shift position secondary_radius = tm.outer_sec_radius(shift_distance) - BELT_HEIGHT / 2 - helix_angle, angle_multiplier, net = _secondary_helix_force_kernel( - torque, - spring_torque_breakdown.net, - secondary_radius, - self.ramp.slope(shift_distance), - ) + # Helix angle at current position + # Negative slope because helix ramps down with increasing shift + helix_angle = np.arctan(-self.ramp.slope(shift_distance)) + + # Convert torque to axial force through helix geometry + # F = (τ + τ_spring) / (2 * tan(α) * r) + angle_multiplier = 2 * np.tan(helix_angle) * secondary_radius + + # Net + net = (torque + spring_torque_breakdown.net) / angle_multiplier return HelixForceBreakdown( feedbackTorque=torque, @@ -285,7 +249,7 @@ def _calculate_spring_tors_torque( """Calculate torsion spring torque (resists helix cam rotation).""" # Clamp shift distance to valid range # TODO: Remove - shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) + shift_distance = np.clip(shift_distance, 0, MAX_SHIFT) # Calculate cam rotation from shift distance (approximation) # TODO: Improve relationship between shift and rotation @@ -313,7 +277,6 @@ def _calculate_rotation(self, shift_distance: float) -> float: Rotation angle [rad] """ # Find x position that corresponds to this height - shift_distance = np.clip(shift_distance, 0, self._max_supported_shift) x_position = self.ramp.find_x_at_height(-shift_distance) # Get slope at that x position return x_position / HELIX_RADIUS diff --git a/cvtModel/src/cvt_simulator/models/slip_model.py b/cvtModel/src/cvt_simulator/models/slip_model.py index 421bf9ea..b7d5bcfb 100644 --- a/cvtModel/src/cvt_simulator/models/slip_model.py +++ b/cvtModel/src/cvt_simulator/models/slip_model.py @@ -1,4 +1,3 @@ -import math import numpy as np from cvt_simulator.models.dataTypes import SlipBreakdown from cvt_simulator.models.external_load_model import LoadModel @@ -20,53 +19,6 @@ from cvt_simulator.utils.numba_utils import maybe_njit -@maybe_njit(cache=True, fastmath=True) -def _relative_speed_kernel( - primary_angular_velocity: float, - secondary_angular_velocity: float, - cvt_ratio: float, -) -> float: - return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) - - -@maybe_njit(cache=True, fastmath=True) -def _coupling_torque_kernel( - relative_speed: float, - torque_demand: float, - t_max_capacity: float, - slip_speed_smoothing: float, -) -> tuple[float, bool]: - coulomb_torque = t_max_capacity * math.tanh(relative_speed / slip_speed_smoothing) - alpha = min(max(abs(relative_speed) / slip_speed_smoothing, 0.0), 1.0) - torque_demand_clamped = min(max(torque_demand, -t_max_capacity), t_max_capacity) - coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque - return coupling_torque, alpha > 0.0 - - -@maybe_njit(cache=True, fastmath=True) -def _torque_demand_kernel( - engine_torque: float, - load_torque: float, - wheel_inertia: float, - wheel_angular_velocity: float, - engine_to_wheel_ratio: float, - engine_to_wheel_ratio_rate_of_change: float, - engine_inertia: float, -) -> float: - eng_term = engine_torque * wheel_inertia - load_term = engine_inertia * load_torque * engine_to_wheel_ratio - shift_term = ( - engine_inertia - * wheel_inertia - * wheel_angular_velocity - * engine_to_wheel_ratio_rate_of_change - ) - - numerator = eng_term + load_term - shift_term - denominator = wheel_inertia + engine_inertia * engine_to_wheel_ratio**2 - return numerator / denominator - - class SlipModel: slip_speed_threshold: float = 2 # rad/s slip_speed_smoothing: float = 5.0 @@ -111,13 +63,24 @@ def get_breakdown(self, state: SystemState) -> SlipBreakdown: state.car_velocity * wheel_to_sec_ratio, tm.current_cvt_ratio(state.shift_distance), ) - coupling_torque, _ = _coupling_torque_kernel( - relative_speed, - torque_demand, - t_max_capacity, - self.slip_speed_smoothing, + + # 1) Smooth Coulomb-like torque based on slip speed + # For large |relative_speed|, tanh -> ±1, so |coupling_torque| -> t_max_capacity + # For small |relative_speed|, torque ~ (t_max_capacity / slip_speed_smoothing) * relative_speed (viscous-ish) + coulomb_torque = t_max_capacity * np.tanh( + relative_speed / self.slip_speed_smoothing ) + # 2) Optionally respect torque_demand near zero slip by blending + # When relative_speed is small, use torque_demand (clamped); + # as slip grows, fade to the Coulomb model. + v_blend = self.slip_speed_smoothing # you can use same scale or a separate one + alpha = np.clip(abs(relative_speed) / v_blend, 0.0, 1.0) + + torque_demand_clamped = np.clip(torque_demand, -t_max_capacity, t_max_capacity) + + coupling_torque = (1.0 - alpha) * torque_demand_clamped + alpha * coulomb_torque + # 3) Define is_slipping for diagnostics (no effect on dynamics) is_slipping = abs(relative_speed) > self.slip_speed_threshold @@ -149,16 +112,21 @@ def get_torque_demand(self, state: SystemState): * GEARBOX_RATIO ) - coupling_torque = _torque_demand_kernel( - engine_torque, - load_torque, - wheel_inertia, - wheel_angular_velocity, - engine_to_wheel_ratio, - engine_to_wheel_ratio_rate_of_change, - ENGINE_INERTIA, + # The secret butter + eng_term = engine_torque * wheel_inertia + load_term = ENGINE_INERTIA * load_torque * engine_to_wheel_ratio + shift_term = ( + ENGINE_INERTIA + * wheel_inertia + * wheel_angular_velocity + * engine_to_wheel_ratio_rate_of_change ) + numerator = eng_term + load_term - shift_term + denominator = wheel_inertia + ENGINE_INERTIA * engine_to_wheel_ratio**2 + + coupling_torque = numerator / denominator + return coupling_torque def get_wheel_speed(self, car_velocity: float): @@ -185,17 +153,14 @@ def calculate_t_max(self, state: SystemState) -> tuple[float, float]: secondary_t_max = max(0, secondary_t_max) return primary_t_max, secondary_t_max + @staticmethod + @maybe_njit(cache=True, fastmath=True) def _relative_speed( - self, primary_angular_velocity: float, secondary_angular_velocity: float, cvt_ratio: float, ) -> float: - return _relative_speed_kernel( - primary_angular_velocity, - secondary_angular_velocity, - cvt_ratio, - ) + return primary_angular_velocity - (secondary_angular_velocity * cvt_ratio) # def _is_slipping( # self,