Propulsion Model Improvements: Efficiency & Autotune Alignment

.gitignore

@@ -11,4 +11,9 @@ venv/
 *_out/
 
 .antigravitycli/
-scratch/
+scratch/
+
+# Test and Coverage
+.coverage
+.pytest_cache/
+htmlcov/

pythrust/propulsion/autotune.py

@@ -17,6 +17,7 @@
 
 from pythrust.propellers.database import PropellerEntry
 from .models import BatterySpec, MotorSpec, PropellerSpec, SystemSpec
+from .solver import evaluate_propulsion_state
 
 
 @dataclass(frozen=True)
@@ -144,18 +145,13 @@ def calibrate(
             )
 
         current_scale = max(p.current_a for p in valid)
-        kt = 60.0 / (2.0 * math.pi * motor.kv_rpm_per_v)
-
         def _motor_state(rpm: float):
-            n = max(rpm / 60.0, 1e-6)
-            j = airspeed_mps / (n * propeller.diameter_m) if propeller.diameter_m > 0 else 0.0
-            ct, cp = prop_entry.get_coefficients(rpm, j)
+            res = evaluate_propulsion_state(motor, propeller, prop_entry, rho, airspeed_mps, rpm)
+            ct, cp, j, torque_nm, i_motor, v_back = res
             if cp <= 0.0 or ct < 0.0 or j < 0.0:
                 return None
-            torque_nm = cp * rho * (n ** 2) * (propeller.diameter_m ** 5) / (2.0 * math.pi)
-            i_motor = torque_nm / kt + motor.no_load_current_a
-            v_back = rpm / motor.kv_rpm_per_v
-            v_m = v_back + i_motor * motor.resistance_ohm
+            v_m = v_back + i_motor * motor.get_winding_resistance(i_motor)
+            n = max(rpm / 60.0, 1e-6)
             thrust_g = ct * rho * (n ** 2) * (propeller.diameter_m ** 4) * 1000.0 / 9.80665
             return i_motor, v_m, thrust_g
 
@@ -168,7 +164,7 @@ def _residuals(params: list) -> list:
                     res.append(1.0)
                     continue
                 i_motor, v_m, _ = state
-                i_bat_pred = (v_m * i_motor + (i_motor ** 2) * r_sys) / battery.voltage_v
+                i_bat_pred = (v_m * i_motor + (i_motor ** 2) * r_sys) / (battery.voltage_v * max(1e-6, battery.discharge_efficiency))
                 res.append((i_bat_pred - pt.current_a) / current_scale)
             return res
 
@@ -194,7 +190,7 @@ def _residuals(params: list) -> list:
                 continue
             i_motor, v_m, thrust_pred_g = state
             thrust_errs.append(thrust_pred_g - pt.thrust_g)
-            i_bat_pred = (v_m * i_motor + (i_motor ** 2) * r_sys_fit) / battery.voltage_v
+            i_bat_pred = (v_m * i_motor + (i_motor ** 2) * r_sys_fit) / (battery.voltage_v * max(1e-6, battery.discharge_efficiency))
             current_errs.append(i_bat_pred - pt.current_a)
 
         rmse_thrust = (

pythrust/propulsion/models.py

@@ -119,3 +119,6 @@ class OperatingPoint:
     motor_voltage_v: float
     is_feasible: bool
     infeasible_reason: Optional[str] = None
+    propeller_efficiency: float = 0.0
+    motor_efficiency: float = 0.0
+    system_efficiency: float = 0.0

pythrust/propulsion/solver.py

@@ -24,6 +24,29 @@ class SolverConfig:
     max_iter: int = 100
 
 
+def evaluate_propulsion_state(
+    motor: MotorSpec,
+    propeller: PropellerSpec,
+    prop_entry: PropellerEntry,
+    rho: float,
+    airspeed_mps: float,
+    rpm: float,
+) -> tuple[float, float, float, float, float, float]:
+    """Calculate aerodynamic and motor electrical states at a given RPM."""
+    n = max(rpm / 60.0, 1e-6)
+    j = airspeed_mps / (n * propeller.diameter_m) if propeller.diameter_m > 0 else 0.0
+    ct, cp = prop_entry.get_coefficients(rpm, j)
+
+    torque_nm = cp * rho * (n**2) * (propeller.diameter_m**5) / (2.0 * math.pi)
+    kt = 30.0 / (math.pi * motor.kv_rpm_per_v * motor.torque_constant_kv_ratio)
+    current_a = torque_nm / kt + motor.get_no_load_current(rpm)
+
+    # Back-EMF with magnetic lag: V_back = (omega * (1 + tau*omega)) / Kv
+    v_back = (rpm / motor.kv_rpm_per_v) * (1.0 + motor.magnetic_lag_tau * rpm * (math.pi / 30.0))
+
+    return ct, cp, j, torque_nm, current_a, v_back
+
+
 class PropulsionSolver:
     """Solve equilibrium RPM for a given operating condition."""
 
@@ -139,18 +162,7 @@ def _evaluate_state(
         airspeed_mps: float,
         rpm: float,
     ) -> tuple[float, float, float, float, float, float]:
-        n = max(rpm / 60.0, 1e-6)
-        j = airspeed_mps / (n * propeller.diameter_m) if propeller.diameter_m > 0 else 0.0
-        ct, cp = prop_entry.get_coefficients(rpm, j)
-
-        torque_nm = cp * rho * (n**2) * (propeller.diameter_m**5) / (2.0 * math.pi)
-        kt = 30.0 / (math.pi * motor.kv_rpm_per_v * motor.torque_constant_kv_ratio)
-        current_a = torque_nm / kt + motor.get_no_load_current(rpm)
-
-        # Back-EMF with magnetic lag: V_back = (omega * (1 + tau*omega)) / Kv
-        v_back = (rpm / motor.kv_rpm_per_v) * (1.0 + motor.magnetic_lag_tau * rpm * (math.pi / 30.0))
-
-        return ct, cp, j, torque_nm, current_a, v_back
+        return evaluate_propulsion_state(motor, propeller, prop_entry, rho, airspeed_mps, rpm)
 
     def _build_point(
         self,
@@ -179,6 +191,22 @@ def _build_point(
         motor_power_w = motor_voltage_v * current_a
         battery_power_w = (motor_power_w + (current_a ** 2) * system.resistance_ohm) / max(1e-6, battery.discharge_efficiency)
 
+        # Efficiency calculations
+        if shaft_power_w > 0.0:
+            propeller_efficiency = (thrust_n * airspeed_mps) / shaft_power_w
+        else:
+            propeller_efficiency = 0.0
+
+        if motor_power_w > 0.0:
+            motor_efficiency = shaft_power_w / motor_power_w
+        else:
+            motor_efficiency = 0.0
+
+        if battery_power_w > 0.0:
+            system_efficiency = (thrust_n * airspeed_mps) / battery_power_w
+        else:
+            system_efficiency = 0.0
+
         feasible = True
         reason = None
         if current_a > motor.current_max_a:
@@ -187,6 +215,11 @@ def _build_point(
         if ct < 0.0 or cp < 0.0 or j < 0.0:
             feasible = False
             reason = "invalid_coefficients"
+        if (propeller_efficiency > 1.0001 or propeller_efficiency < 0.0 or
+            motor_efficiency > 1.0001 or motor_efficiency < 0.0 or
+            system_efficiency > 1.0001 or system_efficiency < 0.0):
+            feasible = False
+            reason = "invalid_efficiency"
 
         return OperatingPoint(
             rpm=float(rpm),
@@ -202,6 +235,9 @@ def _build_point(
             motor_voltage_v=float(motor_voltage_v),
             is_feasible=feasible,
             infeasible_reason=reason,
+            propeller_efficiency=float(propeller_efficiency),
+            motor_efficiency=float(motor_efficiency),
+            system_efficiency=float(system_efficiency),
         )
 
     def _build_infeasible_point(
@@ -240,6 +276,9 @@ def _build_infeasible_point(
             motor_voltage_v=point.motor_voltage_v,
             is_feasible=False,
             infeasible_reason=reason,
+            propeller_efficiency=point.propeller_efficiency,
+            motor_efficiency=point.motor_efficiency,
+            system_efficiency=point.system_efficiency,
         )
 
     @staticmethod

tests/test_models.py

@@ -102,8 +102,34 @@ def test_operating_point():
         battery_power_w=310.0,
         motor_current_a=25.0,
         motor_voltage_v=12.0,
-        is_feasible=True
+        is_feasible=True,
+        propeller_efficiency=0.65,
+        motor_efficiency=0.83,
+        system_efficiency=0.5395
     )
     assert op.rpm == 8000.0
     assert op.is_feasible is True
     assert op.infeasible_reason is None
+    assert op.propeller_efficiency == 0.65
+    assert op.motor_efficiency == 0.83
+    assert op.system_efficiency == 0.5395
+
+    # Test default values
+    op_default = OperatingPoint(
+        rpm=8000.0,
+        advance_ratio=0.4,
+        ct=0.08,
+        cp=0.04,
+        thrust_n=15.0,
+        torque_nm=0.3,
+        shaft_power_w=250.0,
+        motor_power_w=300.0,
+        battery_power_w=310.0,
+        motor_current_a=25.0,
+        motor_voltage_v=12.0,
+        is_feasible=True
+    )
+    assert op_default.propeller_efficiency == 0.0
+    assert op_default.motor_efficiency == 0.0
+    assert op_default.system_efficiency == 0.0
+

tests/test_solver.py

@@ -171,3 +171,95 @@ def test_estimate_j_max(test_setup):
     empty_meta = PropellerMetadata("empty", "A", "B", 10.0, 5.0, 2, "csv")
     empty_entry = PropellerEntry(metadata=empty_meta, data_by_rpm={})
     assert PropulsionSolver._estimate_j_max(empty_entry) == 0.6
+
+
+def test_solver_efficiencies(test_setup):
+    motor, battery, system, propeller, prop_entry = test_setup
+    # Configure higher rpm_min to prevent rpm_min J from exceeding j_max and returning inf
+    solver = PropulsionSolver(SolverConfig(rpm_min=2000.0))
+
+    # 1. Static case (airspeed = 0)
+    op_static = solver.solve_operating_point(
+        motor=motor,
+        battery=battery,
+        system=system,
+        propeller=propeller,
+        prop_entry=prop_entry,
+        rho=1.225,
+        airspeed_mps=0.0,
+        throttle=0.8
+    )
+    assert op_static.is_feasible is True
+    assert op_static.propeller_efficiency == 0.0
+    assert op_static.system_efficiency == 0.0
+    assert 0.0 < op_static.motor_efficiency <= 1.0
+
+    # 2. Dynamic case (airspeed = 5.0 m/s)
+    op_dynamic = solver.solve_operating_point(
+        motor=motor,
+        battery=battery,
+        system=system,
+        propeller=propeller,
+        prop_entry=prop_entry,
+        rho=1.225,
+        airspeed_mps=5.0,
+        throttle=0.8
+    )
+    assert op_dynamic.is_feasible is True
+    # Efficiencies should be positive and physically reasonable
+    assert 0.0 < op_dynamic.propeller_efficiency <= 1.0
+    assert 0.0 < op_dynamic.motor_efficiency <= 1.0
+    assert 0.0 < op_dynamic.system_efficiency <= 1.0
+
+    # System efficiency should be propeller_eff * motor_eff * battery_discharge_eff approximately
+    # Since battery discharge efficiency is 0.98 and sys resistance is small:
+    # eta_sys = (T * V) / P_battery. eta_motor = P_shaft / P_motor. eta_prop = (T * V) / P_shaft.
+    # P_battery = (P_motor + I^2 * R_sys) / eta_batt_discharge.
+    # Let's verify mathematically:
+    expected_sys = op_dynamic.propeller_efficiency * op_dynamic.motor_efficiency
+    # Since there are small electrical line losses (R_sys = 0.015) and battery efficiency (0.98),
+    # actual system efficiency will be slightly lower than expected_sys.
+    assert op_dynamic.system_efficiency <= expected_sys + 1e-5
+
+
+def test_solver_invalid_efficiency(test_setup):
+    motor, battery, system, propeller, _ = test_setup
+
+    # Create custom propeller data where Ct/Cp = 15.0 at J = 0.2
+    # With J = 0.2, eta_prop = (Ct/Cp) * J = 15.0 * 0.2 = 3.0 (> 1.0)
+    meta = PropellerMetadata(
+        id="bad_prop",
+        manufacturer="APC",
+        model="SF",
+        diameter_in=10.0,
+        pitch_in=4.7,
+        blade_count=2,
+        data_csv="bad.csv"
+    )
+    data = {
+        8000.0: [
+            PropellerDataPoint(j=0.0, ct=0.12, cp=0.06),
+            PropellerDataPoint(j=0.2, ct=0.15, cp=0.01),
+            PropellerDataPoint(j=0.8, ct=0.02, cp=0.01)
+        ]
+    }
+    bad_entry = PropellerEntry(metadata=meta, data_by_rpm=data)
+
+    solver = PropulsionSolver(SolverConfig(rpm_min=2000.0))
+    # Solve at airspeed = 5.0 m/s to land J near 0.2
+    op = solver.solve_operating_point(
+        motor=motor,
+        battery=battery,
+        system=system,
+        propeller=propeller,
+        prop_entry=bad_entry,
+        rho=1.225,
+        airspeed_mps=5.0,
+        throttle=0.8
+    )
+
+    # The solver should calculate the point but mark it infeasible
+    assert op.is_feasible is False
+    assert op.infeasible_reason == "invalid_efficiency"
+
+