From 8001d7412062cae90e18668f0e8ce132d56116d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=BCseyin=20Karakaya?= Date: Wed, 3 Jun 2026 18:47:09 +0300 Subject: [PATCH 1/3] feat(propulsion): add propeller, motor, and system efficiency calculations --- .gitignore | 7 ++++- pythrust/propulsion/models.py | 3 +++ pythrust/propulsion/solver.py | 22 +++++++++++++++ tests/test_models.py | 28 +++++++++++++++++++- tests/test_solver.py | 50 +++++++++++++++++++++++++++++++++++ 5 files changed, 108 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index b15cd393..e09fe5c1 100644 --- a/.gitignore +++ b/.gitignore @@ -11,4 +11,9 @@ venv/ *_out/ .antigravitycli/ -scratch/ \ No newline at end of file +scratch/ + +# Test and Coverage +.coverage +.pytest_cache/ +htmlcov/ \ No newline at end of file diff --git a/pythrust/propulsion/models.py b/pythrust/propulsion/models.py index 363810f8..11b9270a 100644 --- a/pythrust/propulsion/models.py +++ b/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 diff --git a/pythrust/propulsion/solver.py b/pythrust/propulsion/solver.py index 841cd53c..03df3245 100644 --- a/pythrust/propulsion/solver.py +++ b/pythrust/propulsion/solver.py @@ -179,6 +179,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: @@ -202,6 +218,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 +259,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 diff --git a/tests/test_models.py b/tests/test_models.py index 78edb496..5680400d 100644 --- a/tests/test_models.py +++ b/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 + diff --git a/tests/test_solver.py b/tests/test_solver.py index 828d2755..b2e94b82 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -171,3 +171,53 @@ 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 + From eafbc4b57b3989b6ab0c10fa4a44298a62d8175d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=BCseyin=20Karakaya?= Date: Wed, 3 Jun 2026 18:47:35 +0300 Subject: [PATCH 2/3] fix(autotune): align motor and battery physics with solver --- pythrust/propulsion/autotune.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pythrust/propulsion/autotune.py b/pythrust/propulsion/autotune.py index f695589b..23c0cffe 100644 --- a/pythrust/propulsion/autotune.py +++ b/pythrust/propulsion/autotune.py @@ -144,7 +144,7 @@ def calibrate( ) current_scale = max(p.current_a for p in valid) - kt = 60.0 / (2.0 * math.pi * motor.kv_rpm_per_v) + kt = 30.0 / (math.pi * motor.kv_rpm_per_v * motor.torque_constant_kv_ratio) def _motor_state(rpm: float): n = max(rpm / 60.0, 1e-6) @@ -153,9 +153,9 @@ def _motor_state(rpm: float): 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 + i_motor = torque_nm / kt + motor.get_no_load_current(rpm) + v_back = (rpm / motor.kv_rpm_per_v) * (1.0 + motor.magnetic_lag_tau * rpm * (math.pi / 30.0)) + v_m = v_back + i_motor * motor.get_winding_resistance(i_motor) thrust_g = ct * rho * (n ** 2) * (propeller.diameter_m ** 4) * 1000.0 / 9.80665 return i_motor, v_m, thrust_g @@ -168,7 +168,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 +194,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 = ( From d188b3296ecfd8f22ddf3b91c32c3df35494d964 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=BCseyin=20Karakaya?= Date: Thu, 4 Jun 2026 01:31:40 +0300 Subject: [PATCH 3/3] feat(solver): add evaluate_propulsion_state function for efficiency calculations fix(tests): update efficiency assertions to allow equality with 1.0 and add invalid efficiency test --- pythrust/propulsion/autotune.py | 12 +++----- pythrust/propulsion/solver.py | 41 +++++++++++++++++++-------- tests/test_solver.py | 50 ++++++++++++++++++++++++++++++--- 3 files changed, 79 insertions(+), 24 deletions(-) diff --git a/pythrust/propulsion/autotune.py b/pythrust/propulsion/autotune.py index 23c0cffe..6ddb1635 100644 --- a/pythrust/propulsion/autotune.py +++ b/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 = 30.0 / (math.pi * motor.kv_rpm_per_v * motor.torque_constant_kv_ratio) - 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.get_no_load_current(rpm) - v_back = (rpm / motor.kv_rpm_per_v) * (1.0 + motor.magnetic_lag_tau * rpm * (math.pi / 30.0)) 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 diff --git a/pythrust/propulsion/solver.py b/pythrust/propulsion/solver.py index 03df3245..9b844c03 100644 --- a/pythrust/propulsion/solver.py +++ b/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, @@ -203,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), diff --git a/tests/test_solver.py b/tests/test_solver.py index b2e94b82..429c73e5 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -192,7 +192,7 @@ def test_solver_efficiencies(test_setup): 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 + assert 0.0 < op_static.motor_efficiency <= 1.0 # 2. Dynamic case (airspeed = 5.0 m/s) op_dynamic = solver.solve_operating_point( @@ -207,9 +207,9 @@ def test_solver_efficiencies(test_setup): ) 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 + 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: @@ -221,3 +221,45 @@ def test_solver_efficiencies(test_setup): # 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" + +