import numpy as np
import numpy.typing as npt
from ...constants import *
from ..vehicle_models import VehicleModel
[docs]
class Motor:
def __init__(
self,
vehicle_model: VehicleModel,
):
self.vehicle_model = vehicle_model
self.params = vehicle_model.pwrtn.motor
self.efficiency_map = None # will contain polyval coefficients
self.rpm_norm_params = None # needed for numerical stability in polyfit
self.torque_norm_params = None # as above
# Precomputed inverse of the DC-power constraint, indexed by rpm.
# Built once at init; rebuild via build_power_limit_table() if motor
# params are mutated at runtime.
self._pow_lim_rpm_grid = None
# _pow_lim_tau_grid[k] is the max motor torque at _pow_lim_rpm_grid[k]
# such that DC bus power τ·ω/η_elec(rpm,τ) does not exceed pow_lim.
self._pow_lim_tau_grid = None
if self.params.use_efficiency_lut:
# Fill out self.efficiency_map
self.fit_efficiency_data(self.params.efficiency_data, 4, 4)
# Note: if the shape is especially weird, may need to increase
# fit degree. Keep even degrees. Keep degree as low as possible
# to limit evaluation complexity
# Invert the DC-power constraint once so the runtime hot path
# is a single np.interp instead of a bisection loop.
self.build_power_limit_table()
[docs]
def get_torque_at_rpm(self, motor_rpm: float) -> float:
"""
Get motor torque at a given RPM using piecewise function.
The torque curve has three regions:
1. Flat region (0 to fw_rpm): constant max torque
2. Linear decay (fw_rpm to max_rpm): torque decreases linearly
3. Zero region (above max_rpm): no torque
Parameters
----------
motor_rpm : float
Motor RPM
Returns
-------
float
Motor torque at the given RPM (Nm)
"""
if motor_rpm < self.params.fw_rpm:
# Pre-field weakening: constant max torque
return self.params.max_torque
elif motor_rpm <= self.params.max_rpm:
# Field weakening: linear decrease from max_torque to fw_torque
# Handle edge case where fw_rpm == max_rpm (no field weakening region)
if self.params.max_rpm == self.params.fw_rpm:
return self.params.max_torque
slope = (self.params.max_torque - self.params.fw_torque) / (
self.params.max_rpm - self.params.fw_rpm
)
return self.params.max_torque - slope * (motor_rpm - self.params.fw_rpm)
else:
# Above max RPM: no torque
return 0.0
[docs]
def calculate_max_ground_force_and_motor_power(self, v: float, ratio: float):
"""
Calculates force and power output by the motor.
The torque the motor can deliver is constrained by two things:
1. the static torque curve τ_curve(rpm), and
2. the DC-bus power limit pow_lim, via ``τ · ω / η_elec(rpm, τ) ≤
pow_lim``, which is circular in τ because η_elec depends on τ.
When use_efficiency_lut is True this is resolved by an inverse
lookup built once at init in build_power_limit_table().
Only chain/diff efficiencies sit between shaft and wheel, so the
wheel torque uses mechanical_efficiency() (not the full chain).
The inverter and motor losses are upstream of the shaft and affect
battery draw, not wheel torque.
Parameters
----------
v : float
Velocity of the vehicle (m/s)
ratio : float
Gear ratio
Returns
-------
ground_force : float
Force at the ground (N)
motor_power : float
Motor shaft power output (W)
"""
motor_rpm = self.v_to_rpm(v, ratio)
omega = rpm_to_rad_s(motor_rpm)
tau_curve = self.get_torque_at_rpm(motor_rpm)
rpm_grid = self._pow_lim_rpm_grid
tau_grid = self._pow_lim_tau_grid
if (
self.params.use_efficiency_lut
and rpm_grid is not None
and tau_grid is not None
):
# DC-bus power constraint, precomputed inverse.
tau_pow = float(np.interp(motor_rpm, rpm_grid, tau_grid))
else:
# Scalar-efficiency fallback: η_elec is constant in τ, so closed form.
eta_elec = self.electrical_efficiency(motor_rpm, 0)
tau_pow = self.params.pow_lim * eta_elec / max(omega, 1e-6)
motor_torque = min(tau_curve, tau_pow)
# Shaft -> wheel: only mechanical components act here.
wheel_torque = motor_torque * ratio * self.mechanical_efficiency()
ground_force = wheel_torque / in_to_m(self.vehicle_model.tires.tire_radius)
return ground_force, motor_torque * omega
[docs]
def v_to_rpm(self, v: float, ratio: float):
"""
Convert vehicle velocity to motor RPM given the gear ratio.
Parameters
----------
v : float
Velocity of the vehicle (m/s)
ratio : float
Gear ratio
Returns
-------
rpm : float
Motor RPM
"""
tire_cir = circumference(
in_to_m(self.vehicle_model.tires.tire_radius)
) # meters
# Calculate index of torque curve given current velocity
rot_per_sec = v / tire_cir
wheel_rpm = rot_per_sec * 60
return wheel_rpm * ratio # rpm_out = rpm_in / ratio
[docs]
def mechanical_efficiency(self):
"""
Loss factor between motor shaft and wheel.
Only the components that physically sit downstream of the shaft act
here: the chain/belt and the diff. The inverter and motor losses are
electrical and do not reduce wheel torque; they show up as extra DC
current draw and are accounted for in electrical_efficiency().
"""
return self.params.chain_efficiency * self.params.diff_efficiency
[docs]
def electrical_efficiency(self, rpm: float, torque: float):
"""
Loss factor between DC bus and motor shaft.
Covers the inverter (moc_efficiency) and the motor itself. When an
efficiency LUT is in use and both rpm and torque are supplied, the
motor's contribution is pulled from the (rpm, torque) map; otherwise
the scalar motor_efficiency is used.
Parameters
----------
rpm : float
Motor RPM (used if LUT is enabled)
torque : float
Motor torque (used if LUT is enabled)
Returns
-------
float
Motor efficiency * MOC efficiency
"""
if self.params.use_efficiency_lut:
eta_motor = self.eval_efficiency_map(rpm, torque)
else:
eta_motor = self.params.motor_efficiency
return self.params.moc_efficiency * eta_motor
[docs]
def powertrain_efficiency(self, rpm: float, torque: float):
"""
Full DC-bus-to-ground efficiency, used by callers converting between
battery power and wheel power (or vice versa).
Equals electrical_efficiency * mechanical_efficiency when in 'indiv'
mode. In 'sys' mode the lumped system_efficiency parameter is used
directly.
Parameters
----------
rpm : float
Motor RPM (used if efficiency_method is 'indiv')
torque : float
Motor torque (used if efficiency_method is 'indiv')
Returns
-------
float
Overall powertrain efficiency (0.0 to 1.0)
"""
# TODO: DIFF EFFICIENCY is set to 100% for now
if self.params.efficiency_method == "sys":
return self.params.system_efficiency
elif self.params.efficiency_method == "indiv":
return (
self.electrical_efficiency(rpm, torque) * self.mechanical_efficiency()
)
[docs]
def build_power_limit_table(self, n_rpm=200):
"""
Precompute the DC-power-limited max motor torque at each rpm.
At a given rpm, the binding DC constraint is
``τ · ω / η_elec(rpm, τ) ≤ pow_lim``,
which is implicit in τ because η_elec depends on τ via the LUT. We
invert this once by bisecting τ in [0, τ_curve(rpm)] at each rpm in
a fixed grid, then look the result up with np.interp at runtime.
At rpm=0 there is no power binding (ω=0 ⇒ P_DC=0 for any τ), so the
static torque curve dominates. Where the static curve already fits
under pow_lim, no inversion is needed.
Rebuild this table if pow_lim, moc_efficiency, or the LUT itself
changes after construction.
Parameters
----------
n_rpm : int
Number of RPM points to precompute between 0 and max_rpm. Default is 200.
"""
rpms = np.linspace(0.0, float(self.params.max_rpm), n_rpm)
tau_pow = np.empty_like(rpms)
pow_lim = self.params.pow_lim
for k, rpm in enumerate(rpms):
tau_curve = self.get_torque_at_rpm(float(rpm))
omega = rpm_to_rad_s(float(rpm))
if omega < 1e-6:
# P_DC = 0 at zero rpm regardless of τ.
tau_pow[k] = tau_curve
continue
# Cheap early-out: if running at full torque is already under
# the power limit, the static curve binds and there's nothing
# to invert.
eta_at_curve = self.electrical_efficiency(float(rpm), tau_curve)
if tau_curve * omega / eta_at_curve <= pow_lim:
tau_pow[k] = tau_curve
continue
# Bisect τ in [0, tau_curve] for τ·ω/η_elec(rpm,τ) = pow_lim.
# P_DC(τ) is non-decreasing in τ because η_elec is bounded
# (eval_efficiency_map clips to [0.66, 0.99]), so the τ term
# dominates as τ → 0+.
lo, hi = 0.0, tau_curve
for _ in range(40): # ample for double precision
mid = 0.5 * (lo + hi)
eta_mid = self.electrical_efficiency(float(rpm), mid)
if mid * omega / eta_mid > pow_lim:
hi = mid
else:
lo = mid
tau_pow[k] = 0.5 * (lo + hi)
self._pow_lim_rpm_grid = rpms
self._pow_lim_tau_grid = tau_pow
[docs]
def fit_efficiency_data(
self, points: list[list[float]], degv: int, degt: int
) -> npt.NDArray:
"""
Fit a 2D polynomial z = f(x, y) to data points [[x, y, z], ...]
Parameters
----------
points : list of lists
Each inner list is [rpm, torque, efficiency] for a data point.
degv : int
Degree of the polynomial in the rpm (x) direction.
degt : int
Degree of the polynomial in the torque (y) direction.
Returns
-------
npt.NDArray
Coefficient matrix C[i,j] for x^i y^j
"""
pts = np.asarray(points, float)
rpm, torque, z = pts.T
# Store normalization parameters
self.rpm_norm_params = (rpm.min(), rpm.max())
self.torque_norm_params = (torque.min(), torque.max())
# Normalize to [0, 1] for numerical stability
x = (rpm - self.rpm_norm_params[0]) / (
self.rpm_norm_params[1] - self.rpm_norm_params[0]
)
y = (torque - self.torque_norm_params[0]) / (
self.torque_norm_params[1] - self.torque_norm_params[0]
)
# Build tensor-product basis
terms = [(x**i) * (y**j) for i in range(degv + 1) for j in range(degt + 1)]
A = np.column_stack(terms)
# Solve least squares
coeffs, *_ = np.linalg.lstsq(A, z, rcond=None)
self.efficiency_map = coeffs.reshape(degv + 1, degt + 1)
return self.efficiency_map
[docs]
def eval_efficiency_map(self, rpm: float, torque: float) -> float:
"""
Evaluate 2D efficiency map (input :math:`(motor rpm, motor torque)`)
with edge clamping and distance penalty for OOB.
Parameters
----------
rpm : float
Motor RPM
torque : float
Motor torque
Returns
-------
float
Motor efficiency (0.66 to 0.99)
"""
if self.efficiency_map is None:
print("WARNING: EFFICIENCY MAP IS EVALUATING BUT DOES NOT EXIST")
return self.params.motor_efficiency
C = self.efficiency_map
rpm_min, rpm_max = self.rpm_norm_params
torque_min, torque_max = self.torque_norm_params
rpm_range = rpm_max - rpm_min
torque_range = torque_max - torque_min
# Distance outside bounds (normalized)
rpm_dist = max(rpm_min - rpm, rpm - rpm_max, 0) / rpm_range
torque_dist = max(torque_min - torque, torque - torque_max, 0) / torque_range
total_dist = (rpm_dist**2 + torque_dist**2) ** 0.5
# Clamp inputs to valid range
rpm_clamped = max(rpm_min, min(rpm_max, rpm))
torque_clamped = max(torque_min, min(torque_max, torque))
# Normalize clamped inputs
x = (rpm_clamped - rpm_min) / rpm_range
y = (torque_clamped - torque_min) / torque_range
# Horner's method evaluation
degx = C.shape[0] - 1
degy = C.shape[1] - 1
tmp = [0.0] * (degx + 1)
for i in range(degx + 1):
eff = C[i, degy]
for j in range(degy - 1, -1, -1):
eff = eff * y + C[i, j]
tmp[i] = eff
eff = tmp[degx]
for i in range(degx - 1, -1, -1):
eff = eff * x + tmp[i]
# Apply distance penalty (exponential decay)
if total_dist > 0:
penalty = 2.71828 ** (-total_dist * 0.3) # ~0.74 at dist=1.0
eff = eff * penalty
# Floor at 66%, cap at 99%
return max(0.66, min(0.99, eff))
[docs]
class Powertrain:
def __init__(
self,
vehicle_model: VehicleModel,
):
self.vehicle_model = vehicle_model
self.params = vehicle_model.pwrtn
# Construct Motor object
self.motor = Motor(vehicle_model=vehicle_model)