diff --git a/README.md b/README.md index d9e4f657..4351e362 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -![PyThrust Banner](https://raw.githubusercontent.com/Setuav/PyThrust/main/docs/images/PyThrust_banner.png) +![PyThrust Banner](docs/images/PyThrust_banner.png) [![CI/CD Pipeline](https://github.com/Setuav/PyThrust/actions/workflows/ci-cd.yml/badge.svg)](https://github.com/Setuav/PyThrust/actions/workflows/ci-cd.yml) [![Docs](https://github.com/Setuav/PyThrust/actions/workflows/docs.yml/badge.svg)](https://github.com/Setuav/PyThrust/actions/workflows/docs.yml) @@ -11,13 +11,13 @@ PyThrust is an open-source framework for electric propulsion system analysis, co-design, and parameter optimization in UAV applications. It can be used for multidisciplinary design optimization (MDO) within OpenMDAO. It includes steady-state performance solvers, auto-tuning calibration tools to fit manufacturer test data, and database search tools to map theoretical designs onto real-world brushless motor and propeller catalogs. -## Design and Analysis Visualization +## Feature Visuals -| 1. Propulsion Co-Design Optimization | 2. Propulsion Calibration & Auto-Tuning | +| System Resistance Calibration | OpenMDAO Hover Co-Design | | :---: | :---: | -| ![Propulsion Co-Design Optimization](https://raw.githubusercontent.com/Setuav/PyThrust/main/docs/images/optimize_and_plot_results.png) | ![Propulsion Calibration & Auto-Tuning Results](https://raw.githubusercontent.com/Setuav/PyThrust/main/docs/images/calibration_results.png) | -| **3. Propeller Aerodynamic Coefficients** | **4. Hover Efficiency Heatmap** | -| ![Propeller Aerodynamic Coefficients](https://raw.githubusercontent.com/Setuav/PyThrust/main/docs/images/propeller_coefficients.png) | ![Hover Efficiency Heatmap](https://raw.githubusercontent.com/Setuav/PyThrust/main/docs/images/efficiency_heatmap.png) | +| ![System Resistance Calibration](docs/images/calibration_results.png) | ![OpenMDAO Hover Co-Design](docs/images/optimize_and_plot_results.png) | +| **Empirical Propeller Database** | **Hover Efficiency Map** | +| ![Empirical Propeller Database](docs/images/propeller_coefficients.png) | ![Hover Efficiency Map](docs/images/efficiency_heatmap.png) | ## Documentation @@ -31,7 +31,7 @@ Key sections: - [Propulsion Solver](https://setuav.github.io/PyThrust/usage/) - [Motor Calibration](https://setuav.github.io/PyThrust/motor_calibration/) - [Examples](https://setuav.github.io/PyThrust/examples/) -- [Mathematical Model](https://setuav.github.io/PyThrust/theory/) +- [Propulsion and Battery Theory](https://setuav.github.io/PyThrust/theory/) - [Component Databases](https://setuav.github.io/PyThrust/databases/) ## License diff --git a/data/batteries/example_liion_cell.json b/data/batteries/example_liion_cell.json new file mode 100644 index 00000000..ed73693a --- /dev/null +++ b/data/batteries/example_liion_cell.json @@ -0,0 +1,15 @@ +{ + "name": "Example Li-ion Cell", + "source": "Synthetic example data for PyThrust tests and examples", + "cell": { + "capacity_ah": 4.2, + "cutoff_voltage_v": 2.5, + "charge_voltage_v": 4.2, + "max_current_a": 20.0 + }, + "curves": { + "dod": [0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0], + "ocv_v": [4.20, 4.08, 3.98, 3.83, 3.70, 3.48, 3.20], + "resistance_ohm": [0.020, 0.021, 0.022, 0.026, 0.031, 0.039, 0.055] + } +} diff --git a/docs/api_reference.md b/docs/api_reference.md index f8ada676..e83b6d84 100644 --- a/docs/api_reference.md +++ b/docs/api_reference.md @@ -18,14 +18,31 @@ This page summarizes the main public classes and helpers. For implementation det Use `get_no_load_current(rpm)` and `get_winding_resistance(current_a)` when evaluating speed-dependent or current-dependent motor behavior. -## Battery, System, and Propeller Specs +## Battery Models + +PyThrust exposes battery models from `pythrust.battery`: + +| Class | Purpose | +|---|---| +| `FixedVoltageBattery` | Historical fixed pack-voltage model with scalar discharge efficiency | +| `RateMapBattery` | Equivalent-circuit model using cell OCV and resistance curves | +| `BatteryState` | State of charge/depth of discharge passed to dynamic battery models | +| `BatteryPoint` | Evaluated terminal voltage, current, power, C-rate, efficiency, and feasibility | + +`BatterySpec` remains available from `pythrust.propulsion` as a compatibility +alias for `FixedVoltageBattery`. New code should import `FixedVoltageBattery` +or `RateMapBattery` directly. + +`RateMapBattery.from_json(path, series=..., parallel=...)` loads one cell +dataset and applies the requested pack topology at runtime. + +## System and Propeller Specs | Class | Purpose | |---|---| -| `BatterySpec` | Pack voltage and discharge efficiency | | `SystemSpec` | Lumped electrical resistance for battery, ESC, wires, and connectors | | `PropellerSpec` | Propeller geometry passed to the solver | -| `OperatingPoint` | Solved RPM, thrust, torque, power, current, voltage, efficiency, and feasibility state | +| `OperatingPoint` | Solved RPM, thrust, torque, motor/battery current, voltage, efficiency, and feasibility state | ## Propulsion Solver @@ -35,6 +52,7 @@ Use `get_no_load_current(rpm)` and `get_winding_resistance(current_a)` when eval point = solver.solve_operating_point( motor=motor, battery=battery, + battery_state=state, # required for RateMapBattery system=system, propeller=propeller, prop_entry=prop_entry, @@ -44,6 +62,9 @@ point = solver.solve_operating_point( ) ``` +`battery_state` may be omitted for `FixedVoltageBattery`. It is required for +dynamic battery models such as `RateMapBattery`. + `SolverConfig` controls numerical behavior: | Field | Default | Description | @@ -54,6 +75,17 @@ point = solver.solve_operating_point( | `eps_v` | `1e-8` | Voltage residual tolerance | | `max_iter` | `100` | Maximum root-finder iterations | +`OperatingPoint` includes propulsion outputs such as `rpm`, `thrust_n`, +`torque_nm`, `motor_current_a`, and `motor_voltage_v`, plus battery outputs: + +| Field | Description | +|---|---| +| `battery_power_w` | Battery-side power draw | +| `battery_voltage_v` | Battery terminal pack voltage | +| `battery_current_a` | Battery DC current draw | +| `battery_c_rate` | Cell C-rate for rate-map batteries, or `0.0` for fixed-voltage batteries | +| `battery_efficiency` | Battery model discharge efficiency at the solved point | + ## Propeller Database `PropellerDatabase` loads JSON metadata and CSV performance tables: @@ -129,4 +161,4 @@ system = result.to_system_spec() `pythrust.openmdao.PropulsionComponent` wraps `PropulsionSolver` as an `ExplicitComponent` for optimization models. -Inputs include motor parameters, battery voltage, system resistance, propeller diameter, throttle, density, and airspeed. Outputs include RPM, thrust, torque, battery current, battery power, motor current, motor voltage, and feasibility. +Inputs include motor parameters, fixed battery voltage, system resistance, propeller diameter, throttle, density, and airspeed. Outputs include RPM, thrust, torque, battery current, battery power, motor current, motor voltage, and feasibility. diff --git a/docs/battery_model.md b/docs/battery_model.md new file mode 100644 index 00000000..1a6ca3cd --- /dev/null +++ b/docs/battery_model.md @@ -0,0 +1,268 @@ +# Battery Model + +Status: implemented model note for issue #3. + +PyThrust supports a fixed-voltage battery model and a lightweight rate-map +battery model. The fixed-voltage path is useful for quick propulsion sizing, +but it hides two effects that matter for electric aircraft performance studies: + +- The terminal voltage drops with load. +- The usable energy depends on discharge rate and state of charge. + +This page defines the rate-map model now implemented in PyThrust. The model is +inspired by Robert A. McDonald's +`bat-perf` model and the paper "Battery Knockdown Factors for Conceptual +Design". + +!!! abstract "Model scope" + PyThrust uses a compact equivalent-circuit battery model for sizing and optimization. It is intended to capture voltage sag and usable energy trends without introducing electrochemical simulation inputs. + +## Goals + +The model is intended to: + +- Stay fast enough for sizing sweeps, optimizers, and OpenMDAO workflows. +- Use manufacturer-accessible data such as capacity, voltage limits, current + limits, discharge curves, and C-rate maps. +- Support point-performance analysis at a specified state of charge. +- Support mission integration by advancing battery state through time. +- Preserve a simple fixed-voltage battery path for examples and low-fidelity + studies. + +The model does not try to replace electrochemical simulation tools. PyThrust +does not need electrode-level parameters, diffusion constants, thermal cell +models, or PyBaMM-style electrochemistry for this feature. + +## Core Equations + +The paper treats battery state with two governing equations: conservation of +charge and an algebraic terminal-voltage equation. + +Depth of discharge is: + +$$ +x = 1 - z +$$ + +where $z$ is state of charge. With discharge current $I$ taken as positive, +charge conservation is: + +$$ +\frac{dx}{dt} = \frac{I}{Q} +$$ + +where $Q$ is rated cell capacity in ampere-seconds. + +The static equivalent circuit is an open-circuit voltage source in series with +an internal resistance: + +$$ +V(x, I) = OCV(x) - R(x)I +$$ + +where: + +- $V$ is terminal cell voltage. +- $OCV(x)$ is open-circuit cell voltage as a function of depth of discharge. +- $R(x)$ is effective cell resistance as a function of depth of discharge. +- $I$ is cell current, positive during discharge. + +For a pack with cells in series and parallel: + +$$ +V_{pack} = N_s V_{cell} +$$ + +$$ +I_{cell} = \frac{I_{pack}}{N_p} +$$ + +$$ +Q_{pack} = N_p Q_{cell} +$$ + +## bat-perf Mapping + +The `bat-perf` MATLAB model uses the same core variables: + +| PyThrust concept | bat-perf name | Meaning | +|---|---|---| +| `dod` | `dod` | Depth of discharge, from 0 to 1 | +| `capacity_as` | `Qmax` | Cell capacity in ampere-seconds | +| `rated_current_a` | `irated` | 1C current | +| `ocv(dod)` | `OCVfun(dod)` | Open-circuit cell voltage | +| `resistance(dod)` | `Rssfun(dod)` | Steady-state internal resistance | +| `cutoff_voltage_v` | `Vcutoff` | Minimum terminal cell voltage | +| `charge_voltage_v` | `Vcharge` | Maximum terminal cell voltage | + +The most important point-state functions are: + +| Mode | bat-perf function | Equation | +|---|---|---| +| Specified current | `cellStateI` | $V = OCV - RI$ | +| Specified C-rate | `cellStateC` | $I = C Q_{Ah}$, then $V = OCV - RI$ | +| Specified voltage | `cellStateV` | $I = (OCV - V) / R$ | +| Specified power | `cellStateP` | solve $P = I(OCV - RI)$ | +| Specified load resistance | `cellStateR` | $I = OCV / (R + R_{load})$ | + +For specified power, the current is obtained from the quadratic: + +$$ +RI^2 - OCV I + P = 0 +$$ + +using the lower-current discharge root: + +$$ +I = \frac{OCV - \sqrt{OCV^2 - 4RP}}{2R} +$$ + +The maximum deliverable power at a given DOD occurs when the discriminant is +zero: + +$$ +P_{max}(x) = \frac{OCV(x)^2}{4R(x)} +$$ + +The implementation reports infeasible states when requested power exceeds this +limit, when current exceeds the configured limit, or when terminal voltage falls +below cutoff. + +## Python API + +Use explicit names for the two battery fidelities: + +=== "Fixed voltage" + + ```python + from pythrust.battery import FixedVoltageBattery + + battery = FixedVoltageBattery(voltage_v=14.8) + ``` + +=== "Rate map" + + ```python + from pythrust.battery import BatteryState, RateMapBattery + + state = BatteryState(soc=1.0) + battery = RateMapBattery.from_json( + "data/batteries/example_liion_cell.json", + series=4, + parallel=2, + ) + ``` + +!!! tip "Choosing a battery model" + Use `FixedVoltageBattery` for early propulsion sizing and simple examples. Use `RateMapBattery` when load-dependent voltage, C-rate limits, or state of charge affect the result. + +`BatterySpec` is too general once multiple battery models exist. It remains as +a compatibility alias: + +```python +BatterySpec = FixedVoltageBattery +``` + +The alias keeps old code working during the transition, but new code and +documentation should use `FixedVoltageBattery`. + +The common battery behavior is intentionally small: + +```python +voltage = battery.terminal_voltage(current_a=current, state=state) +power = battery.terminal_power(current_a=current, state=state) +``` + +`RateMapBattery` also supports state advancement: + +```python +next_state = battery.step_power(power_w=power, dt_s=dt, state=state) +``` + +For the fixed-voltage model, `terminal_voltage` returns the configured voltage. +For the rate-map model, it evaluates the cell/pack equivalent circuit. + +## Dataset Shape + +The JSON dataset describes one cell. Pack topology is passed when loading the +dataset: + +```json +{ + "name": "Example Li-ion Cell", + "source": "Synthetic example data for PyThrust tests and examples", + "cell": { + "capacity_ah": 4.2, + "cutoff_voltage_v": 2.5, + "charge_voltage_v": 4.2, + "max_current_a": 20.0 + }, + "curves": { + "dod": [0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0], + "ocv_v": [4.20, 4.08, 3.98, 3.83, 3.70, 3.48, 3.20], + "resistance_ohm": [0.020, 0.021, 0.022, 0.026, 0.031, 0.039, 0.055] + } +} +``` + +```python +battery = RateMapBattery.from_json(cell_path, series=4, parallel=2) +``` + +The current implementation interpolates `OCV(dod)` and `R(dod)` directly. A +later calibration utility can derive these curves from manufacturer C-rate +discharge maps. Manufacturer discharge curves are usually terminal voltage +under load, so real datasets should document how `OCV(dod)` and `R(dod)` were +derived. + +## Solver Integration + +The propulsion solver starts from the PWM voltage relation: + +$$ +V_{applied} = throttle \times V_{pack} +$$ + +For a rate-map battery, `V_pack` depends on current and state: + +$$ +V_{applied} = throttle \times V_{pack}(x, I) +$$ + +The root equation therefore becomes: + +$$ +F(RPM) = +V_{back}(RPM) ++ I(RPM)R_{motor} ++ I(RPM)R_{system} +- throttle \times V_{pack}(x, I(RPM)) +$$ + +This keeps the propeller/motor equilibrium as a one-dimensional root solve +because current remains a function of RPM through the propeller torque demand. + +For mission simulation, evaluate each time step with the current state, compute +pack current/power, then advance DOD: + +$$ +x_{next} = x + \frac{I_{cell}}{Q_{cell}} \Delta t +$$ + +## Implementation Status + +The initial implementation includes: + +- `pythrust.battery.FixedVoltageBattery` +- `pythrust.battery.RateMapBattery` +- `pythrust.battery.BatteryState` and `BatteryPoint` +- JSON cell datasets with explicit series and parallel counts at load time +- Solver integration through `solve_operating_point(..., battery_state=...)` +- `OperatingPoint` battery outputs for voltage, current, C-rate, and efficiency +- A runnable rate-map mission example + +## References + +- Robert A. McDonald, "Battery Knockdown Factors for Conceptual Design", + AIAA Aviation Forum, 2024, DOI: `10.2514/6.2024-3903`. +- `ramcdona/bat-perf`: diff --git a/docs/development.md b/docs/development.md index 6e16133c..1261621d 100644 --- a/docs/development.md +++ b/docs/development.md @@ -22,15 +22,24 @@ The CI workflow runs the test suite on Python 3.10, 3.11, and 3.12. ## Run Example Workflows ```bash -PYTHONPATH=. python examples/select_motor.py -PYTHONPATH=. python examples/calibrate_from_datasheet.py -PYTHONPATH=. python examples/optimize_and_plot_propulsion.py +PYTHONPATH=. python examples/calibrate_system_resistance.py +PYTHONPATH=. python examples/rate_map_battery_point_states.py +PYTHONPATH=. python examples/rate_map_battery_mission.py +PYTHONPATH=. python examples/select_motor_from_database.py +PYTHONPATH=. python examples/openmdao_hover_optimization.py ``` See [Examples](examples.md) for a user-facing walkthrough of each script. Generated plots are written under `docs/images/` and are used by the documentation site. +Regenerate the README and documentation feature visuals: + +```bash +PYTHONPATH=. python scripts/generate_readme_plots.py +PYTHONPATH=. python examples/openmdao_hover_optimization.py +``` + ## Documentation Site PyThrust uses MkDocs Material for a clean, searchable static documentation site. diff --git a/docs/examples.md b/docs/examples.md index 1be010e5..c3d9d9a0 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -1,22 +1,25 @@ # Examples -PyThrust includes runnable examples that show the main workflows: solving against catalog data, calibrating losses from measurements, and using OpenMDAO for propulsion co-design. +PyThrust includes runnable examples that show the main workflows: solving and +calibrating propulsion models, inspecting rate-map battery behavior, and using +OpenMDAO for hover co-design. -Run examples from the repository root so relative `data/` and `docs/images/` paths resolve correctly. +Run examples from the repository root so relative `data/` and `docs/images/` +paths resolve correctly. ```bash PYTHONPATH=. python examples/.py ``` ---- - ## Requirements | Example | Extra dependencies | |---|---| -| `calibrate_from_datasheet.py` | Core PyThrust dependencies | -| `select_motor.py` | `openmdao` | -| `optimize_and_plot_propulsion.py` | `openmdao`, `matplotlib` | +| `calibrate_system_resistance.py` | Core PyThrust dependencies | +| `rate_map_battery_point_states.py` | Core PyThrust dependencies | +| `rate_map_battery_mission.py` | Core PyThrust dependencies | +| `select_motor_from_database.py` | `openmdao` | +| `openmdao_hover_optimization.py` | `openmdao`, `matplotlib` | Install the full example environment: @@ -24,17 +27,16 @@ Install the full example environment: pip install -e .[plot,openmdao] ``` ---- - -## Datasheet Calibration +## System Resistance Calibration Script: ```bash -PYTHONPATH=. python examples/calibrate_from_datasheet.py +PYTHONPATH=. python examples/calibrate_system_resistance.py ``` -This example identifies the lumped system resistance for a motor, propeller, battery, ESC, and wiring setup. +This example identifies the lumped system resistance for a motor, propeller, +fixed-voltage battery, ESC, and wiring setup. It uses: @@ -42,7 +44,7 @@ It uses: |---|---| | Motor | Datasheet Kv, resistance, no-load current, and current limit | | Propeller | `APC_13x6.5E` from `data/propellers/apc_202602` | -| Battery | 4S nominal voltage, `14.8 V` | +| Battery | Fixed 4S nominal voltage, `14.8 V` | | Test table | RPM, thrust in grams, and battery current in amps | The output reports: @@ -55,18 +57,66 @@ The output reports: | Thrust R2 | Fit quality for the aerodynamic thrust prediction | | Per-point table | Predicted vs measured thrust/current for each RPM row | -See [Motor Calibration](motor_calibration.md) for the calibration model and equations. +See [Motor Calibration](motor_calibration.md) for the calibration model and +equations. + +![System resistance calibration](images/calibration_results.png) + +## Rate-Map Battery Point States + +Script: + +```bash +PYTHONPATH=. python examples/rate_map_battery_point_states.py +``` -![Calibration results](images/calibration_results.png) +This example loads `data/batteries/example_liion_cell.json`, applies a `4S2P` +pack topology, and evaluates the same state of charge under several requests. ---- +It demonstrates: + +| Query | Meaning | +|---|---| +| `state_at_current` | Terminal voltage and power at a requested pack current | +| `state_at_c_rate` | Current/voltage behavior at a requested cell C-rate | +| `state_at_voltage` | Current required to hold a requested pack voltage | +| `state_at_power` | Current and voltage for a requested pack power | +| `state_at_load_resistance` | Battery behavior under a resistive load | +| Infeasible power | How the model reports a power limit | + +This example is intentionally independent from the propulsion solver. It +validates the battery model surface before mission or solver integration. + +## Rate-Map Battery Mission + +Script: + +```bash +PYTHONPATH=. python examples/rate_map_battery_mission.py +``` + +This example couples `RateMapBattery` to `PropulsionSolver` over a short +segment schedule. Each segment solves the propulsion operating point using the +current battery state, reports pack current and voltage, then advances state of +charge from the solved battery current. + +It demonstrates: + +| Step | Meaning | +|---|---| +| Load cell data | Use `data/batteries/example_liion_cell.json` with explicit series and parallel counts | +| Solve segment | Pass `battery_state` into `solve_operating_point(...)` | +| Read outputs | Inspect `battery_voltage_v` and `battery_current_a` on `OperatingPoint` | +| Advance state | Use `step_current(...)` to update SoC for the next segment | + +![Rate-map battery mission simulation](images/rate_map_battery_mission.png) ## Motor Selection Script: ```bash -PYTHONPATH=. python examples/select_motor.py +PYTHONPATH=. python examples/select_motor_from_database.py ``` This example combines theoretical co-design with real motor database lookup. @@ -74,36 +124,27 @@ This example combines theoretical co-design with real motor database lookup. Workflow: 1. Load `APC_13x6.5E` propeller data. -2. Use OpenMDAO to find an efficient theoretical motor/propeller/throttle combination for hover. +2. Use OpenMDAO to find an efficient theoretical motor/propeller/throttle + combination for hover. 3. Load the brushless motor database from `data/motors`. 4. Search real motors near the optimized Kv and current requirement. 5. Print the top candidates sorted by winding resistance and weight. The optimization target is a hover thrust of `4.903 N`, approximately `500 gf`. -Typical output includes: - -| Output | Meaning | -|---|---| -| Target Kv | Ideal speed constant from the theoretical optimization | -| Target diameter | Optimized propeller diameter | -| Target hover current | Current at the optimized hover point | -| Minimum hover power | Battery power objective value | -| Top motor matches | Closest catalog motors with Kv, resistance, weight, and current limit | - -See [Component Databases](databases.md) for motor catalog format and query helpers. - ---- +See [Component Databases](databases.md) for motor catalog format and query +helpers. -## Propulsion Optimization and Plotting +## OpenMDAO Hover Optimization Script: ```bash -PYTHONPATH=. python examples/optimize_and_plot_propulsion.py +PYTHONPATH=. python examples/openmdao_hover_optimization.py ``` -This example demonstrates OpenMDAO-based propulsion co-design and a parametric Kv sweep. +This example demonstrates OpenMDAO-based propulsion co-design and a parametric +Kv sweep. It performs three stages: @@ -121,9 +162,10 @@ The plot shows: | Panel | Shows | |---|---| -| Power and propeller size vs motor Kv | Hover battery power and optimized propeller diameter | -| Throttle and RPM vs motor Kv | Optimized throttle setting and shaft speed | +| Power and propeller sizing | Hover battery power and optimized propeller diameter | +| Control setting and shaft speed | Optimized throttle setting and shaft speed | -![Propulsion co-design optimization](images/optimize_and_plot_results.png) +![OpenMDAO hover co-design](images/optimize_and_plot_results.png) -See [Propulsion Solver](usage.md) for the operating-point solver used inside the OpenMDAO component. +See [Propulsion Solver](usage.md) for the operating-point solver used inside +the OpenMDAO component. diff --git a/docs/getting_started.md b/docs/getting_started.md index 2e88f4c7..0aa6a208 100644 --- a/docs/getting_started.md +++ b/docs/getting_started.md @@ -48,9 +48,9 @@ The core workflow is: ```python from pathlib import Path +from pythrust.battery import FixedVoltageBattery from pythrust.propellers import PropellerDatabase from pythrust.propulsion import ( - BatterySpec, MotorSpec, PropellerSpec, PropulsionSolver, @@ -67,7 +67,7 @@ motor = MotorSpec( no_load_current_a=1.3, current_max_a=65.0, ) -battery = BatterySpec(voltage_v=14.8) +battery = FixedVoltageBattery(voltage_v=14.8) system = SystemSpec(resistance_ohm=0.05) propeller = PropellerSpec(diameter_m=0.3302, blade_count=2) @@ -92,9 +92,11 @@ print(point.is_feasible) ## Run the Examples ```bash -PYTHONPATH=. python examples/select_motor.py -PYTHONPATH=. python examples/calibrate_from_datasheet.py -PYTHONPATH=. python examples/optimize_and_plot_propulsion.py +PYTHONPATH=. python examples/calibrate_system_resistance.py +PYTHONPATH=. python examples/rate_map_battery_point_states.py +PYTHONPATH=. python examples/rate_map_battery_mission.py +PYTHONPATH=. python examples/select_motor_from_database.py +PYTHONPATH=. python examples/openmdao_hover_optimization.py ``` The plotting and OpenMDAO examples require the optional dependencies shown above. diff --git a/docs/images/calibration_results.png b/docs/images/calibration_results.png index 89284238..bdc925f1 100644 Binary files a/docs/images/calibration_results.png and b/docs/images/calibration_results.png differ diff --git a/docs/images/efficiency_heatmap.png b/docs/images/efficiency_heatmap.png index 89b7461a..1215b6e7 100644 Binary files a/docs/images/efficiency_heatmap.png and b/docs/images/efficiency_heatmap.png differ diff --git a/docs/images/optimize_and_plot_results.png b/docs/images/optimize_and_plot_results.png index be0f074d..ff504873 100644 Binary files a/docs/images/optimize_and_plot_results.png and b/docs/images/optimize_and_plot_results.png differ diff --git a/docs/images/propeller_coefficients.png b/docs/images/propeller_coefficients.png index 9a098ad5..04549ffa 100644 Binary files a/docs/images/propeller_coefficients.png and b/docs/images/propeller_coefficients.png differ diff --git a/docs/images/pybamm_mission_results.png b/docs/images/pybamm_mission_results.png deleted file mode 100644 index 991f5394..00000000 Binary files a/docs/images/pybamm_mission_results.png and /dev/null differ diff --git a/docs/images/rate_map_battery_mission.png b/docs/images/rate_map_battery_mission.png new file mode 100644 index 00000000..4c5797ed Binary files /dev/null and b/docs/images/rate_map_battery_mission.png differ diff --git a/docs/index.md b/docs/index.md index 0a15cf6a..e68804d9 100644 --- a/docs/index.md +++ b/docs/index.md @@ -8,33 +8,78 @@ PyThrust combines empirical propeller data, brushless motor models, battery/syst --- -## Why PyThrust? +## Core Capabilities Electric UAV propulsion design usually crosses several domains: aerodynamics, motor electrical behavior, battery loading, component catalogs, and mission constraints. PyThrust keeps those pieces in one workflow: -* **Coupled operating-point solver:** Solve equilibrium RPM, thrust, torque, current, power, and efficiency for a motor, propeller, battery, and flight condition. -* **Catalog-backed selection:** Query real brushless motor and propeller datasets instead of optimizing against abstract components only. -* **Calibration tools:** Fit lumped system resistance from test-stand data to account for ESC, battery, wiring, and connector losses. -* **Optimization-ready components:** Use the propulsion solver inside OpenMDAO for multidisciplinary design optimization studies. -* **Optimization-ready examples:** Generate plots for design sweeps, calibration quality, propeller coefficients, and hover efficiency maps. +
+ +- **Operating-point solver** + + Solve equilibrium RPM, thrust, torque, current, voltage, power, and efficiency for a coupled motor-propeller-battery system. + +- **Catalog-backed selection** + + Query empirical propeller data and brushless motor records instead of sizing only against abstract component assumptions. + +- **Battery-aware analysis** + + Use fixed-voltage batteries for quick studies or rate-map batteries when voltage sag and state of charge matter. + +- **Calibration and optimization** + + Fit lumped system resistance from test data and use the solver inside OpenMDAO co-design workflows. + +
--- -## Design and Analysis Visualization +## Feature Visuals -| Propulsion Co-Design Optimization | Propulsion Calibration & Auto-Tuning | +| System Resistance Calibration | OpenMDAO Hover Co-Design | | :---: | :---: | -| ![Propulsion co-design optimization](images/optimize_and_plot_results.png) | ![Propulsion calibration and auto-tuning results](images/calibration_results.png) | -| **Propeller Aerodynamic Coefficients** | **Hover Efficiency Heatmap** | -| ![Propeller aerodynamic coefficients](images/propeller_coefficients.png) | ![Hover efficiency heatmap](images/efficiency_heatmap.png) | - -## Key Documentation Sections - -* [**Getting Started**](getting_started.md): Installation, optional extras, and a first operating-point solve. -* [**Propulsion Solver**](usage.md): Solver configuration, feasibility rules, and usage examples. -* [**Motor Calibration**](motor_calibration.md): Fit system resistance from manufacturer or thrust-stand data. -* [**Examples**](examples.md): Runnable scripts for calibration, motor selection, and OpenMDAO optimization. -* [**API Reference**](api_reference.md): Main classes, database loaders, calibration objects, and OpenMDAO wrapper. -* [**Mathematical Model**](theory.md): Propeller, motor, electrical loss, and coupled equilibrium equations. -* [**Component Databases**](databases.md): Propeller CSV/JSON and motor catalog formats. -* [**Development & Testing**](development.md): Local setup, tests, examples, and documentation build commands. +| ![System resistance calibration](images/calibration_results.png) | ![OpenMDAO hover co-design](images/optimize_and_plot_results.png) | +| **Empirical Propeller Database** | **Hover Efficiency Map** | +| ![Empirical propeller database](images/propeller_coefficients.png) | ![Hover efficiency map](images/efficiency_heatmap.png) | + +## Explore the Docs + +
+ +- **Getting Started** + + Install PyThrust, add optional extras, and run the first operating-point solve. + + [Open guide](getting_started.md) + +- **Propulsion Solver** + + Review solver inputs, equations, feasibility rules, result fields, and usage examples. + + [Open guide](usage.md) + +- **Battery Model** + + Understand fixed-voltage and rate-map batteries, JSON datasets, and solver integration. + + [Open reference](battery_model.md) + +- **Examples** + + Run calibration, motor selection, battery, and OpenMDAO workflows from the repository. + + [Open examples](examples.md) + +- **Theory** + + Trace the propeller, motor, electrical loss, and battery equations used by the solver. + + [Open theory](theory.md) + +- **API Reference** + + Find the main classes, database loaders, calibration objects, and OpenMDAO wrapper. + + [Open API](api_reference.md) + +
diff --git a/docs/motor_calibration.md b/docs/motor_calibration.md index 6d66da5c..ae1f6503 100644 --- a/docs/motor_calibration.md +++ b/docs/motor_calibration.md @@ -181,11 +181,12 @@ The calibration outcome returns residuals and $R^2$ values to evaluate propeller ## Usage ```python +from pythrust.battery import FixedVoltageBattery from pythrust.propulsion.autotune import ManufacturerTestPoint, PropulsionCalibrator -from pythrust.propulsion import MotorSpec, BatterySpec, SystemSpec, PropellerSpec +from pythrust.propulsion import MotorSpec, SystemSpec, PropellerSpec motor = MotorSpec(kv_rpm_per_v=860, resistance_ohm=0.0258, no_load_current_a=1.3, current_max_a=65) -battery = BatterySpec(voltage_v=14.8) +battery = FixedVoltageBattery(voltage_v=14.8) system = SystemSpec(resistance_ohm=0.05) # starting guess propeller = PropellerSpec(diameter_m=0.3302) diff --git a/docs/theory.md b/docs/theory.md index f82e39f2..1c76a889 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -1,4 +1,4 @@ -# Propulsion System Mathematical Model +# Propulsion and Battery Model Theory This document defines the core mathematical model for the electric propulsion system (propeller + brushless motor + battery/ESC). @@ -23,6 +23,9 @@ This document defines the core mathematical model for the electric propulsion sy | $R_{\text{system}}$ | Lumped transmission system resistance | $\Omega$ | | $V_m$ | Motor terminal voltage | $\text{V}$ | | $V_{\text{back}}$ | Back-EMF voltage | $\text{V}$ | +| $V_{\text{pack}}$ | Battery terminal pack voltage | $\text{V}$ | +| $x$ | Battery depth of discharge | Dimensionless | +| $N_s, N_p$ | Pack cells in series and parallel | Dimensionless | --- @@ -74,7 +77,39 @@ where $\tau$ is the magnetic lag time constant (set to $0$ for first-order model --- -## 4) System Electrical Resistance & Power Chain +## 4) Battery Model + +The fixed-voltage battery model uses a constant terminal pack voltage: + +$$ +V_{\text{pack}} = V_{\text{nominal}} +$$ + +The rate-map battery model uses a cell open-circuit voltage curve and a cell +resistance curve indexed by depth of discharge: + +$$ +V_{\text{cell}}(x, I_{\text{cell}}) = +OCV(x) - R(x)I_{\text{cell}} +$$ + +For a pack topology: + +$$ +V_{\text{pack}} = N_s V_{\text{cell}} +$$ + +$$ +I_{\text{cell}} = \frac{I_{\text{pack}}}{N_p} +$$ + +During mission stepping, state advances by charge conservation: + +$$ +x_{\text{next}} = x + \frac{I_{\text{cell}}}{Q_{\text{cell}}}\Delta t +$$ + +## 5) System Electrical Resistance & Power Chain The voltage drops across ESC MOSFETs, battery internal resistance, cables, and connectors are modeled as a lumped transmission system resistance ($R_{\text{system}}$): @@ -88,14 +123,23 @@ $$ P_{\text{battery}} = \frac{V_m I + I^2 R_{\text{system}}}{\eta_{\text{discharge}}} $$ +For `RateMapBattery`, the equivalent-circuit terminal voltage already captures +the cell-level voltage sag from the battery curves. `R_system` should still be +used for wiring, ESC, connector, or other losses outside the battery model. + --- -## 5) Coupled Equilibrium Condition +## 6) Coupled Equilibrium Condition For a given throttle setting and airspeed, the equilibrium shaft speed (RPM) is the solution of the coupled electrical and aerodynamic torque equilibrium equation: $$ -F(\text{RPM}) = \text{throttle} \times V_{\text{pack}} - \Big( V_{\text{back}}(\text{RPM}) + I(\text{RPM}) (R + R_{\text{system}}) \Big) = 0 +F(\text{RPM}) = +V_{\text{back}}(\text{RPM}) ++ I(\text{RPM})R ++ I(\text{RPM})R_{\text{system}} +- \text{throttle} \times V_{\text{pack}}(x, I(\text{RPM})) += 0 $$ A root-finding method (e.g., Brent's method) solves $F(\text{RPM}) = 0$ for RPM. Once the equilibrium RPM is determined, $T, Q, P_{\text{shaft}}, I$, and efficiency parameters are calculated. @@ -111,3 +155,7 @@ A root-finding method (e.g., Brent's method) solves $F(\text{RPM}) = 0$ for RPM. 2. **Second-Order DC Electric Motor Model** Mark Drela, MIT Aero & Astro, March 2006 [PDF Link](https://web.mit.edu/drela/Public/web/qprop/motor2_theory.pdf) + +3. **Battery Knockdown Factors for Conceptual Design** + Robert A. McDonald, AIAA Aviation Forum, 2024 + DOI: `10.2514/6.2024-3903` diff --git a/docs/usage.md b/docs/usage.md index 7e7383f1..d0dfb4fd 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -2,9 +2,11 @@ This guide describes how PyThrust solves a motor-propeller operating point and how to interpret the result. -The solver answers one core question: +!!! question "Core solver question" + Given a motor, propeller, battery, throttle, and airspeed, what shaft RPM makes the electrical motor model and propeller aerodynamic load agree? -> Given a motor, propeller, battery, throttle, and airspeed, what shaft RPM makes the electrical motor model and propeller aerodynamic load agree? +!!! abstract "Use this page when" + You are wiring the solver into an analysis script, checking why a point is infeasible, or confirming how battery voltage enters the RPM root solve. --- @@ -15,7 +17,7 @@ The operating-point solver combines five inputs: | Input | Class or value | Purpose | |---|---|---| | Motor | `MotorSpec` | Kv, winding resistance, no-load current, current limit, optional higher-order loss terms | -| Battery | `BatterySpec` | Pack voltage and discharge efficiency | +| Battery | `FixedVoltageBattery` or `RateMapBattery` | Fixed pack voltage or state-dependent pack voltage from cell curves | | System | `SystemSpec` | Lumped resistance for ESC, battery internal resistance, wiring, and connectors | | Propeller | `PropellerSpec` + `PropellerEntry` | Geometry plus empirical aerodynamic coefficients | | Flight condition | `rho`, `airspeed_mps`, `throttle` | Air density, freestream speed, and commanded voltage fraction | @@ -43,11 +45,31 @@ The operating-point solver combines five inputs: | $V_{\text{back}}$ | Motor back-EMF voltage | | $V_m$ | Motor terminal voltage | | $V_{\text{pack}}$ | Battery pack voltage | +| $x$ | Battery depth of discharge | --- ## What the Solver Does +```mermaid +flowchart LR + motor["MotorSpec"]:::motorStyle + prop["PropellerSpec + PropellerEntry"]:::propStyle + battery["FixedVoltageBattery or RateMapBattery"]:::batteryStyle + system["SystemSpec"] + flight["rho, airspeed, throttle"] + solver["PropulsionSolver"] + point["OperatingPoint"] + + motor --> solver + prop --> solver + battery --> solver + system --> solver + flight --> solver + solver --> point + point --> outputs["RPM, thrust, current, voltage, efficiency"] +``` + ### Step 1: Evaluate the propeller at a candidate RPM For each candidate RPM, PyThrust converts shaft speed to revolutions per second: @@ -129,7 +151,7 @@ $$ | `resistance_quadratic` | Adds current-dependent winding resistance | | `iron_loss_exponent` | Scales no-load current with RPM using a power law | -When these fields are left at their defaults, the model reduces to the simpler first-order motor equations above. See [Mathematical Model](theory.md) for the full equations and Drela references. +When these fields are left at their defaults, the model reduces to the simpler first-order motor equations above. See [Propulsion and Battery Theory](theory.md) for the full equations and Drela references. --- @@ -142,6 +164,17 @@ V_{\text{applied}} = \text{throttle} \cdot V_{\text{pack}} $$ +For `FixedVoltageBattery`, `V_pack` is the configured pack voltage. For +`RateMapBattery`, `V_pack` is evaluated from the current battery state and the +current implied by the candidate RPM: + +!!! note "Rate-map batteries are state dependent" + A `RateMapBattery` solve needs a `BatteryState` because terminal voltage depends on both depth of discharge and current. + +$$ +V_{\text{pack}} = V_{\text{pack}}(x, I) +$$ + The solver searches for the RPM where the motor voltage demand plus system voltage drop equals the applied voltage: $$ @@ -168,7 +201,7 @@ Before solving, PyThrust builds a search interval: | Bound | How it is chosen | |---|---| | Lower RPM | At least `SolverConfig.rpm_min`; raised if airspeed and propeller `J` range require a higher RPM | -| Upper RPM | Estimated from `motor.kv_rpm_per_v * battery.voltage_v * throttle`, then expanded by `rpm_max_margin` | +| Upper RPM | Estimated from `motor.kv_rpm_per_v * battery.terminal_voltage(0 A) * throttle`, then expanded by `rpm_max_margin` | If the residual does not change sign inside the bracket, the point is returned as infeasible with reason `no_bracket`. @@ -202,6 +235,10 @@ The numerical behavior of the root finder is controlled by `SolverConfig`: | `shaft_power_w` | Mechanical shaft power | | `motor_power_w` | Electrical power at motor terminals | | `battery_power_w` | Battery-side power including system losses | +| `battery_voltage_v` | Battery terminal pack voltage at the solved point | +| `battery_current_a` | Battery DC current draw at the solved point | +| `battery_c_rate` | Cell C-rate for rate-map batteries, or `0.0` for fixed-voltage batteries | +| `battery_efficiency` | Discharge efficiency from the active battery model | | `motor_current_a` | Motor winding current | | `motor_voltage_v` | Motor terminal voltage | | `propeller_efficiency` | Propulsive efficiency based on thrust power | @@ -222,6 +259,10 @@ An operating point is marked as infeasible when one of these checks fails: | `no_bracket` | The RPM search interval does not contain a valid voltage-balance root | | `no_convergence` | The root finder did not converge | | `current_limit` | `motor_current_a` exceeds `current_max_a` | +| `battery_current_limit` | Rate-map battery current exceeds the configured cell current limit | +| `battery_voltage_cutoff` | Rate-map battery terminal cell voltage falls below cutoff | +| `battery_voltage_limit` | Rate-map battery terminal cell voltage exceeds the configured charge voltage | +| `battery_state_limit` | Rate-map battery state is outside the valid range | | `invalid_coefficients` | Propeller coefficient lookup produced non-physical values | | `invalid_efficiency` | Computed efficiency is outside the physically expected range | @@ -234,9 +275,9 @@ Here is a complete example showing how to load a propeller dataset, define speci ```python from pathlib import Path +from pythrust.battery import FixedVoltageBattery from pythrust.propellers import PropellerDatabase from pythrust.propulsion import ( - BatterySpec, MotorSpec, PropellerSpec, PropulsionSolver, @@ -254,7 +295,7 @@ motor = MotorSpec( current_max_a=65.0, ) -battery = BatterySpec( +battery = FixedVoltageBattery( voltage_v=14.8, discharge_efficiency=1.0, ) @@ -277,7 +318,35 @@ point = solver.solve_operating_point( print(f"RPM : {point.rpm:.1f}") print(f"Thrust : {point.thrust_n:.2f} N") print(f"Motor Current : {point.motor_current_a:.2f} A") +print(f"Battery Voltage : {point.battery_voltage_v:.2f} V") +print(f"Battery Current : {point.battery_current_a:.2f} A") print(f"Battery Power : {point.battery_power_w:.1f} W") print(f"Feasible : {point.is_feasible}") print(f"Reason : {point.infeasible_reason}") ``` + +For a rate-map battery, load a cell dataset, create an explicit battery state, +and pass that state into the solver: + +```python +from pythrust.battery import BatteryState, RateMapBattery + +battery = RateMapBattery.from_json( + "data/batteries/example_liion_cell.json", + series=4, + parallel=2, +) +state = BatteryState(soc=0.85) + +point = solver.solve_operating_point( + motor=motor, + battery=battery, + battery_state=state, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=10.0, + throttle=0.6, +) +``` diff --git a/examples/calibrate_from_datasheet.py b/examples/calibrate_system_resistance.py similarity index 62% rename from examples/calibrate_from_datasheet.py rename to examples/calibrate_system_resistance.py index 4c6a2822..1d9f8977 100644 --- a/examples/calibrate_from_datasheet.py +++ b/examples/calibrate_system_resistance.py @@ -1,47 +1,52 @@ -"""Example: calibrate system resistance from a manufacturer thrust/current table. +"""Calibrate system resistance from a manufacturer thrust/current table. This example demonstrates how to use PropulsionCalibrator to identify the -system resistance parameter (lumped ESC, cable, and battery resistance) so that -PyThrust's model matches the manufacturer's performance data. +system resistance parameter so PyThrust's fixed-voltage propulsion model +matches the manufacturer's performance data. Usage:: - python examples/calibrate_from_datasheet.py + PYTHONPATH=. python examples/calibrate_system_resistance.py """ import math +import sys from pathlib import Path +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +from pythrust.battery import FixedVoltageBattery from pythrust.propellers import PropellerDatabase from pythrust.propulsion import ( - BatterySpec, MotorSpec, PropellerSpec, SystemSpec, ) from pythrust.propulsion.autotune import ManufacturerTestPoint, PropulsionCalibrator -# ── 1. Motor parameters — taken directly from the manufacturer datasheet ────── +# 1. Motor parameters taken directly from the manufacturer datasheet. motor = MotorSpec( kv_rpm_per_v=860.0, # RPM/V - resistance_ohm=0.0258, # Ω (terminal resistance) - no_load_current_a=1.3, # A (idle current at rated voltage) - current_max_a=65.0, # A (continuous current limit) + resistance_ohm=0.0258, # Ohm, terminal resistance + no_load_current_a=1.3, # A, idle current at rated voltage + current_max_a=65.0, # A, continuous current limit ) -# ── 2. Battery and propeller specs ──────────────────────────────────────────── -battery = BatterySpec(voltage_v=14.8) # 4S LiPo at nominal voltage -system = SystemSpec(resistance_ohm=0.05) # initial guess (will be fitted) -propeller = PropellerSpec(diameter_m=0.3302) # 13-inch propeller +# 2. Fixed-voltage battery, initial system resistance, and propeller geometry. +battery = FixedVoltageBattery(voltage_v=14.8) # 4S LiPo at nominal voltage +system = SystemSpec(resistance_ohm=0.05) # initial guess, fitted below +propeller = PropellerSpec(diameter_m=0.3302) # 13-inch propeller -# ── 3. Load propeller aerodynamic data ──────────────────────────────────────── +# 3. Load propeller aerodynamic data. db = PropellerDatabase() db.load(Path("data/propellers/apc_202602"), strict=False) prop_entry = db.get("APC_13x6.5E") if prop_entry is None: raise SystemExit("Propeller 'APC_13x6.5E' not found in dataset.") -# ── 4. Manufacturer test table (RPM, Thrust in grams, Battery Current in Amps) ─ +# 4. Manufacturer test table: RPM, thrust in grams, and battery current in amps. RAW_TABLE = [ {"rpm": 3897, "thrust_g": 500, "current_a": 3.9}, {"rpm": 4804, "thrust_g": 750, "current_a": 6.7}, @@ -64,30 +69,30 @@ for row in RAW_TABLE ] -# ── 5. Calibrate ────────────────────────────────────────────────────────────── +# 5. Calibrate the lumped system resistance. calibrator = PropulsionCalibrator() result = calibrator.calibrate( points, motor, battery, system, propeller, prop_entry, rho=1.225 ) -print("─── Calibration Result ──────────────────────────────────────") +print("--- Calibration Result --------------------------------------") print(f" System resistance : {result.system_resistance_ohm:.5f} ohm") print(f" Thrust RMSE : {result.rmse_thrust_g:.1f} g") print(f" Current RMSE : {result.rmse_current_a:.2f} A") -print(f" R² (thrust) : {result.r_squared_thrust:.4f}") +print(f" R^2 (thrust) : {result.r_squared_thrust:.4f}") print(f" Converged : {result.converged}") print(f" Points used : {result.n_points}") for w in result.warnings: - print(f" ⚠ {w}") + print(f" WARNING: {w}") print() -# ── 6. Validate: compare predictions against the table at each RPM ──────────── +# 6. Validate by comparing predictions against the table at each RPM. system_calibrated = result.to_system_spec() kt = 60.0 / (2.0 * math.pi * motor.kv_rpm_per_v) header = f"{'RPM':>6} {'Thrust Pred':>12} {'Thrust Meas':>12} {'Thrust Err':>11} {'Current Pred':>13} {'Current Meas':>13} {'Current Err':>11}" print(header) -print("─" * len(header)) +print("-" * len(header)) for pt in points: n = pt.rpm / 60.0 ct, cp = prop_entry.get_coefficients(pt.rpm, 0.0) diff --git a/examples/openmdao_hover_optimization.py b/examples/openmdao_hover_optimization.py new file mode 100644 index 00000000..e77dbdc6 --- /dev/null +++ b/examples/openmdao_hover_optimization.py @@ -0,0 +1,219 @@ +"""OpenMDAO hover co-design optimization and Kv sweep. + +This example optimizes motor Kv, propeller diameter, and throttle for a fixed +hover thrust, then sweeps Kv to show the design trend. + +Usage:: + + PYTHONPATH=. python examples/openmdao_hover_optimization.py +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np +import openmdao.api as om + +from pythrust.openmdao import PropulsionComponent +from pythrust.propellers import PropellerDatabase + + +COLORS = { + "blue": "#2563eb", + "green": "#059669", + "orange": "#d97706", + "purple": "#7c3aed", + "red": "#dc2626", + "muted": "#6b7280", +} + +FIGSIZE = (9.6, 5.4) +DPI = 170 + + +def main(): + # 1. Load propeller aerodynamic data. + db = PropellerDatabase() + db.load(Path("data/propellers/apc_202602"), strict=False) + prop_entry = db.get("APC_13x6.5E") + if prop_entry is None: + print("Error: Propeller 'APC_13x6.5E' not found.") + sys.exit(1) + + # 2. Set up the baseline optimization problem. + prob = om.Problem() + model = prob.model + + comp = PropulsionComponent(prop_entry=prop_entry) + model.add_subsystem("prop", comp, promotes=["*"]) + + model.add_design_var("kv_rpm_per_v", lower=600.0, upper=1200.0, ref=1000.0) + model.add_design_var("diameter_m", lower=0.254, upper=0.381, ref=0.3) + model.add_design_var("throttle", lower=0.2, upper=0.9, ref=0.5) + + model.add_constraint("thrust_n", equals=4.903, ref=5.0) + model.add_constraint("motor_current_a", upper=25.0, ref=25.0) + model.add_objective("battery_power_w", ref=100.0) + + prob.driver = om.ScipyOptimizeDriver(optimizer="SLSQP") + prob.driver.options["disp"] = False + prob.setup() + + prob.set_val("kv_rpm_per_v", 860.0) + prob.set_val("resistance_ohm", 0.0258) + prob.set_val("no_load_current_a", 1.3) + prob.set_val("current_max_a", 65.0) + prob.set_val("voltage_v", 14.8) + prob.set_val("resistance_system_ohm", 0.095) + prob.set_val("diameter_m", 0.3302) + prob.set_val("throttle", 0.5) + prob.set_val("rho", 1.225) + prob.set_val("airspeed_mps", 0.0) + + print("=== Running Baseline Model ===") + prob.run_model() + base_kv = float(prob.get_val("kv_rpm_per_v")[0]) + base_dia = float(prob.get_val("diameter_m")[0]) * 39.3701 + base_power = float(prob.get_val("battery_power_w")[0]) + base_throttle = float(prob.get_val("throttle")[0]) * 100.0 + base_rpm = float(prob.get_val("rpm")[0]) + print(f" Motor Kv : {base_kv:.1f} RPM/V") + print(f" Propeller Diameter: {base_dia:.2f} inches") + print(f" Hover Throttle : {base_throttle:.1f}%") + print(f" Hover Shaft Speed : {base_rpm:.0f} RPM") + print(f" Battery Power Draw: {base_power:.2f} W") + + print("\n=== Running Co-Design Optimization ===") + prob.run_driver() + opt_kv = float(prob.get_val("kv_rpm_per_v")[0]) + opt_dia = float(prob.get_val("diameter_m")[0]) * 39.3701 + opt_power = float(prob.get_val("battery_power_w")[0]) + opt_throttle = float(prob.get_val("throttle")[0]) * 100.0 + opt_rpm = float(prob.get_val("rpm")[0]) + print(f" Optimal Motor Kv : {opt_kv:.1f} RPM/V") + print(f" Optimal Diameter : {opt_dia:.2f} inches") + print(f" Optimal Throttle : {opt_throttle:.1f}%") + print(f" Optimal Shaft Speed : {opt_rpm:.0f} RPM") + print(f" Minimum Power Draw : {opt_power:.2f} W") + print(f" Power Reduction : {base_power - opt_power:.2f} W ({(base_power - opt_power)/base_power*100:.1f}%)") + + # 3. Sweep Kv and re-optimize diameter/throttle at each point. + print("\n=== Performing Kv Parametric Sweep ===") + kv_sweep = np.linspace(600, 1200, 13) + power_draws = [] + optimal_diameters = [] + optimal_throttles = [] + optimal_rpms = [] + + sweep_prob = om.Problem() + sweep_model = sweep_prob.model + sweep_model.add_subsystem("prop", PropulsionComponent(prop_entry=prop_entry), promotes=["*"]) + sweep_model.add_design_var("diameter_m", lower=0.254, upper=0.381, ref=0.3) + sweep_model.add_design_var("throttle", lower=0.2, upper=0.9, ref=0.5) + sweep_model.add_constraint("thrust_n", equals=4.903, ref=5.0) + sweep_model.add_objective("battery_power_w", ref=100.0) + sweep_prob.driver = om.ScipyOptimizeDriver(optimizer="SLSQP") + sweep_prob.driver.options["disp"] = False + sweep_prob.setup() + + sweep_prob.set_val("resistance_ohm", 0.0258) + sweep_prob.set_val("no_load_current_a", 1.3) + sweep_prob.set_val("current_max_a", 65.0) + sweep_prob.set_val("voltage_v", 14.8) + sweep_prob.set_val("resistance_system_ohm", 0.095) + sweep_prob.set_val("rho", 1.225) + sweep_prob.set_val("airspeed_mps", 0.0) + + for temp_kv in kv_sweep: + sweep_prob.set_val("kv_rpm_per_v", temp_kv) + sweep_prob.set_val("diameter_m", 0.3302) + sweep_prob.set_val("throttle", 0.5) + sweep_prob.run_driver() + + power = float(sweep_prob.get_val("battery_power_w")[0]) + dia = float(sweep_prob.get_val("diameter_m")[0]) * 39.3701 + thr = float(sweep_prob.get_val("throttle")[0]) * 100.0 + rpm = float(sweep_prob.get_val("rpm")[0]) + + power_draws.append(power) + optimal_diameters.append(dia) + optimal_throttles.append(thr) + optimal_rpms.append(rpm) + + # 4. Plot optimization trends. + plt.rcParams.update( + { + "figure.facecolor": "white", + "axes.facecolor": "white", + "axes.edgecolor": "#374151", + "grid.color": "#d1d5db", + "font.size": 10, + "axes.titlesize": 12, + "figure.titlesize": 15, + } + ) + fig, (ax1, ax3) = plt.subplots(1, 2, figsize=FIGSIZE) + + # Left plot: Power and Diameter vs Motor Kv + ax1.grid(True) + ax1.set_xlabel("Motor Kv [RPM/V]") + ax1.set_ylabel("Hover battery power [W]", color=COLORS["blue"]) + ax1.plot(kv_sweep, power_draws, color=COLORS["blue"], marker="o", linewidth=2.1, label="power") + ax1.tick_params(axis="y", labelcolor=COLORS["blue"]) + + # Highlight the global optimum found in step 2 + ax1.plot(opt_kv, opt_power, marker="*", color=COLORS["red"], markersize=13, label="best point") + + ax2 = ax1.twinx() + ax2.set_ylabel("Optimized propeller diameter [in]", color=COLORS["orange"]) + ax2.plot(kv_sweep, optimal_diameters, color=COLORS["orange"], marker="s", linestyle="--", linewidth=2.0, label="diameter") + ax2.tick_params(axis="y", labelcolor=COLORS["orange"]) + + ax1.set_title("Power and propeller sizing") + + # Right plot: Throttle and RPM vs Motor Kv + ax3.grid(True) + ax3.set_xlabel("Motor Kv [RPM/V]") + ax3.set_ylabel("Hover throttle [%]", color=COLORS["green"]) + ax3.plot(kv_sweep, optimal_throttles, color=COLORS["green"], marker="^", linewidth=2.1, label="throttle") + ax3.tick_params(axis="y", labelcolor=COLORS["green"]) + + ax4 = ax3.twinx() + ax4.set_ylabel("Hover shaft speed [RPM]", color=COLORS["purple"]) + ax4.plot(kv_sweep, optimal_rpms, color=COLORS["purple"], marker="d", linestyle="--", linewidth=2.0, label="RPM") + ax4.tick_params(axis="y", labelcolor=COLORS["purple"]) + + ax3.set_title("Control setting and shaft speed") + + fig.suptitle("OpenMDAO Hover Co-Design", fontweight="bold") + fig.text( + 0.5, + 0.92, + "Optimize motor Kv, propeller diameter, and throttle for 500 g hover", + ha="center", + color=COLORS["muted"], + fontsize=10, + ) + fig.subplots_adjust( + left=0.09, + right=0.90, + bottom=0.17, + top=0.78, + wspace=0.48, + ) + + output_image = Path("docs/images/optimize_and_plot_results.png") + plt.savefig(output_image, dpi=DPI) + print(f"\nSaved default style plot to: {output_image.resolve()}") + +if __name__ == "__main__": + main() diff --git a/examples/optimize_and_plot_propulsion.py b/examples/optimize_and_plot_propulsion.py deleted file mode 100644 index 575f6717..00000000 --- a/examples/optimize_and_plot_propulsion.py +++ /dev/null @@ -1,181 +0,0 @@ -"""OpenMDAO Propulsion Co-Design Optimization and Sweep. - -This script demonstrates how to use OpenMDAO to optimize a drone propulsion system -(Motor Kv, Propeller Diameter, and Throttle) for maximum efficiency in hover, and -plots the design space sweep using the default Matplotlib appearance. -""" - -import sys -from pathlib import Path -import numpy as np -import matplotlib.pyplot as plt -import openmdao.api as om - -# Append project root to sys.path -sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) - -from pythrust.propellers import PropellerDatabase -from pythrust.openmdao import PropulsionComponent - - -def main(): - # 1. Load Propeller Database - db = PropellerDatabase() - db.load(Path("data/propellers/apc_202602"), strict=False) - prop_entry = db.get("APC_13x6.5E") - if prop_entry is None: - print("Error: Propeller 'APC_13x6.5E' not found.") - sys.exit(1) - - # 2. Setup the Optimization Problem - prob = om.Problem() - model = prob.model - - # Add the PyThrust Propulsion component - comp = PropulsionComponent(prop_entry=prop_entry) - model.add_subsystem('prop', comp, promotes=['*']) - - # Define Design Variables (with scaling ref) - model.add_design_var('kv_rpm_per_v', lower=600.0, upper=1200.0, ref=1000.0) - model.add_design_var('diameter_m', lower=0.254, upper=0.381, ref=0.3) - model.add_design_var('throttle', lower=0.2, upper=0.9, ref=0.5) - - # Define Constraints and Objective - model.add_constraint('thrust_n', equals=4.903, ref=5.0) # 500 gf thrust target - model.add_constraint('motor_current_a', upper=25.0, ref=25.0) - model.add_objective('battery_power_w', ref=100.0) - - # Setup Scipy SLSQP Driver - prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP') - prob.driver.options['disp'] = False # quiet run - prob.setup() - - # Set initial baseline conditions - prob.set_val('kv_rpm_per_v', 860.0) - prob.set_val('resistance_ohm', 0.0258) - prob.set_val('no_load_current_a', 1.3) - prob.set_val('current_max_a', 65.0) - prob.set_val('voltage_v', 14.8) - prob.set_val('resistance_system_ohm', 0.095) - prob.set_val('diameter_m', 0.3302) - prob.set_val('throttle', 0.5) - prob.set_val('rho', 1.225) - prob.set_val('airspeed_mps', 0.0) - - print("=== Running Baseline Model ===") - prob.run_model() - base_kv = float(prob.get_val('kv_rpm_per_v')[0]) - base_dia = float(prob.get_val('diameter_m')[0]) * 39.3701 - base_power = float(prob.get_val('battery_power_w')[0]) - base_throttle = float(prob.get_val('throttle')[0]) * 100.0 - base_rpm = float(prob.get_val('rpm')[0]) - print(f" Motor Kv : {base_kv:.1f} RPM/V") - print(f" Propeller Diameter: {base_dia:.2f} inches") - print(f" Hover Throttle : {base_throttle:.1f}%") - print(f" Hover Shaft Speed : {base_rpm:.0f} RPM") - print(f" Battery Power Draw: {base_power:.2f} W") - - print("\n=== Running Co-Design Optimization ===") - prob.run_driver() - opt_kv = float(prob.get_val('kv_rpm_per_v')[0]) - opt_dia = float(prob.get_val('diameter_m')[0]) * 39.3701 - opt_power = float(prob.get_val('battery_power_w')[0]) - opt_throttle = float(prob.get_val('throttle')[0]) * 100.0 - opt_rpm = float(prob.get_val('rpm')[0]) - print(f" Optimal Motor Kv : {opt_kv:.1f} RPM/V") - print(f" Optimal Diameter : {opt_dia:.2f} inches") - print(f" Optimal Throttle : {opt_throttle:.1f}%") - print(f" Optimal Shaft Speed : {opt_rpm:.0f} RPM") - print(f" Minimum Power Draw : {opt_power:.2f} W") - print(f" Power Reduction : {base_power - opt_power:.2f} W ({(base_power - opt_power)/base_power*100:.1f}%)") - - # 3. Perform Parametric Sweep of Kv - print("\n=== Performing Kv Parametric Sweep ===") - kv_sweep = np.linspace(600, 1200, 13) - power_draws = [] - optimal_diameters = [] - optimal_throttles = [] - optimal_rpms = [] - - # Setup a sub-problem for the sweeps - sweep_prob = om.Problem() - sweep_model = sweep_prob.model - sweep_model.add_subsystem('prop', PropulsionComponent(prop_entry=prop_entry), promotes=['*']) - sweep_model.add_design_var('diameter_m', lower=0.254, upper=0.381, ref=0.3) - sweep_model.add_design_var('throttle', lower=0.2, upper=0.9, ref=0.5) - sweep_model.add_constraint('thrust_n', equals=4.903, ref=5.0) - sweep_model.add_objective('battery_power_w', ref=100.0) - sweep_prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP') - sweep_prob.driver.options['disp'] = False - sweep_prob.setup() - - sweep_prob.set_val('resistance_ohm', 0.0258) - sweep_prob.set_val('no_load_current_a', 1.3) - sweep_prob.set_val('current_max_a', 65.0) - sweep_prob.set_val('voltage_v', 14.8) - sweep_prob.set_val('resistance_system_ohm', 0.095) - sweep_prob.set_val('rho', 1.225) - sweep_prob.set_val('airspeed_mps', 0.0) - - for temp_kv in kv_sweep: - sweep_prob.set_val('kv_rpm_per_v', temp_kv) - sweep_prob.set_val('diameter_m', 0.3302) - sweep_prob.set_val('throttle', 0.5) - sweep_prob.run_driver() - - power = float(sweep_prob.get_val('battery_power_w')[0]) - dia = float(sweep_prob.get_val('diameter_m')[0]) * 39.3701 - thr = float(sweep_prob.get_val('throttle')[0]) * 100.0 - rpm = float(sweep_prob.get_val('rpm')[0]) - - power_draws.append(power) - optimal_diameters.append(dia) - optimal_throttles.append(thr) - optimal_rpms.append(rpm) - - # 4. Plot using Default Matplotlib Appearance - fig, (ax1, ax3) = plt.subplots(1, 2, figsize=(12, 5)) - - # Left plot: Power and Diameter vs Motor Kv - ax1.grid(True) - ax1.set_xlabel('Motor Kv (RPM/V)') - ax1.set_ylabel('Hover Battery Power (W)', color='C0') - ax1.plot(kv_sweep, power_draws, color='C0', marker='o', label='Power (W)') - ax1.tick_params(axis='y', labelcolor='C0') - - # Highlight the global optimum found in step 2 - ax1.plot(opt_kv, opt_power, marker='*', color='red', markersize=12, label='Global Optimum') - - ax2 = ax1.twinx() - ax2.set_ylabel('Optimal Propeller Diameter (in)', color='C1') - ax2.plot(kv_sweep, optimal_diameters, color='C1', marker='s', linestyle='--', label='Diameter (in)') - ax2.tick_params(axis='y', labelcolor='C1') - - ax1.set_title('Power and Propeller Size vs Motor Kv') - - # Right plot: Throttle and RPM vs Motor Kv - ax3.grid(True) - ax3.set_xlabel('Motor Kv (RPM/V)') - ax3.set_ylabel('Hover Throttle (%)', color='C2') - ax3.plot(kv_sweep, optimal_throttles, color='C2', marker='^', label='Throttle (%)') - ax3.tick_params(axis='y', labelcolor='C2') - - ax4 = ax3.twinx() - ax4.set_ylabel('Hover RPM', color='C4') - ax4.plot(kv_sweep, optimal_rpms, color='C4', marker='d', linestyle='--', label='RPM') - ax4.tick_params(axis='y', labelcolor='C4') - - ax3.set_title('Throttle and RPM vs Motor Kv') - - fig.suptitle('Propulsion Co-Design Optimization', fontsize=14) - plt.tight_layout() - - output_image = Path("docs/images/optimize_and_plot_results.png") - plt.savefig(output_image, bbox_inches='tight') - print(f"\nSaved default style plot to: {output_image.resolve()}") - - plt.show() - - -if __name__ == '__main__': - main() diff --git a/examples/rate_map_battery_mission.py b/examples/rate_map_battery_mission.py new file mode 100644 index 00000000..20a47acd --- /dev/null +++ b/examples/rate_map_battery_mission.py @@ -0,0 +1,110 @@ +"""Simulate a simple mission with a rate-map battery. + +This example couples RateMapBattery to PropulsionSolver. Each segment solves +the propulsion operating point from the current battery state, then advances +state of charge using the solved battery current. + +Usage:: + + PYTHONPATH=. python examples/rate_map_battery_mission.py +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +from pythrust.battery import BatteryState, RateMapBattery +from pythrust.propellers import PropellerDatabase +from pythrust.propulsion import MotorSpec, PropellerSpec, PropulsionSolver, SystemSpec + + +MISSION_SEGMENTS = [ + {"name": "takeoff", "duration_s": 45.0, "throttle": 0.70, "airspeed_mps": 0.0}, + {"name": "climb", "duration_s": 90.0, "throttle": 0.65, "airspeed_mps": 5.0}, + {"name": "cruise", "duration_s": 180.0, "throttle": 0.50, "airspeed_mps": 10.0}, + {"name": "return", "duration_s": 120.0, "throttle": 0.45, "airspeed_mps": 8.0}, +] + + +def load_propeller(): + db = PropellerDatabase() + db.load(Path("data/propellers/apc_202602"), strict=False) + prop_entry = db.get("APC_13x6.5E") + if prop_entry is None: + raise SystemExit("Propeller 'APC_13x6.5E' not found in dataset.") + return prop_entry + + +def main(): + battery_dataset = Path("data/batteries/example_liion_cell.json") + battery = RateMapBattery.from_json(battery_dataset, series=4, parallel=2) + state = BatteryState(soc=0.95) + + motor = MotorSpec( + kv_rpm_per_v=860.0, + resistance_ohm=0.0258, + no_load_current_a=1.3, + current_max_a=65.0, + ) + system = SystemSpec(resistance_ohm=0.095) + propeller = PropellerSpec(diameter_m=0.3302) + prop_entry = load_propeller() + solver = PropulsionSolver() + + total_energy_wh = 0.0 + total_time_s = 0.0 + + print("Rate-map battery mission") + print(f"Pack: {battery.name}, {battery.series}S{battery.parallel}P") + print(f"Initial SoC: {state.soc * 100.0:.1f}%") + print() + + header = f"{'Segment':<10}{'min':>6}{'thr':>7}{'SoC':>17}{'I [A]':>9}{'V [V]':>9}" + print(header) + print("-" * len(header)) + + for segment in MISSION_SEGMENTS: + op = solver.solve_operating_point( + motor=motor, + battery=battery, + battery_state=state, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=segment["airspeed_mps"], + throttle=segment["throttle"], + ) + + if not op.is_feasible: + reason = op.infeasible_reason or "unknown" + raise SystemExit(f"Mission segment '{segment['name']}' is infeasible: {reason}") + + next_state = battery.step_current( + state=state, + current_a=op.battery_current_a, + dt_s=segment["duration_s"], + ) + + total_energy_wh += op.battery_power_w * segment["duration_s"] / 3600.0 + total_time_s += segment["duration_s"] + + print( + f"{segment['name']:<10}" + f"{segment['duration_s'] / 60.0:>6.1f}" + f"{segment['throttle'] * 100.0:>6.0f}%" + f"{state.soc * 100.0:>7.1f}% -> {next_state.soc * 100.0:>5.1f}%" + f"{op.battery_current_a:>9.2f}" + f"{op.battery_voltage_v:>9.2f}" + ) + state = next_state + + print() + print(f"Summary: {total_time_s / 60.0:.1f} min, {total_energy_wh:.1f} Wh, final SoC {state.soc * 100.0:.1f}%") + + +if __name__ == "__main__": + main() diff --git a/examples/rate_map_battery_point_states.py b/examples/rate_map_battery_point_states.py new file mode 100644 index 00000000..711a9669 --- /dev/null +++ b/examples/rate_map_battery_point_states.py @@ -0,0 +1,62 @@ +"""Evaluate rate-map battery point states. + +This example loads a cell-level rate-map dataset, applies a 4S2P pack topology, +and evaluates the same battery state under current, C-rate, voltage, power, and +load-resistance requests. + +Usage:: + + PYTHONPATH=. python examples/rate_map_battery_point_states.py +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +from pythrust.battery import BatteryState, RateMapBattery + + +def print_point(label, point): + """Print one compact battery point-state row.""" + status = "OK" if point.is_feasible else f"NO ({point.infeasible_reason})" + print( + f"{label:<18}" + f"{point.terminal_voltage_v:>10.2f}" + f"{point.current_a:>10.2f}" + f"{point.power_w:>11.2f}" + f"{point.c_rate:>9.2f}" + f"{point.efficiency:>9.3f}" + f" {status}" + ) + + +def main(): + dataset = Path("data/batteries/example_liion_cell.json") + battery = RateMapBattery.from_json(dataset, series=4, parallel=2) + state = BatteryState(soc=0.60) + + print("Rate-map battery point states") + print(f"Cell dataset : {dataset}") + print(f"Pack topology: {battery.series}S{battery.parallel}P") + print(f"State : SoC={state.soc:.2f}, DOD={state.dod:.2f}") + print() + print(f"{'Case':<18}{'Vpack':>10}{'Ipack':>10}{'Ppack':>11}{'C-rate':>9}{'eta':>9} Status") + print("-" * 79) + + print_point("current 12 A", battery.state_at_current(state=state, current_a=12.0)) + print_point("1.5 C", battery.state_at_c_rate(state=state, c_rate=1.5)) + print_point("voltage 14 V", battery.state_at_voltage(state=state, voltage_v=14.0)) + print_point("power 180 W", battery.state_at_power(state=state, power_w=180.0)) + print_point("load 1.5 ohm", battery.state_at_load_resistance(state=state, resistance_ohm=1.5)) + print_point("too much power", battery.state_at_power(state=state, power_w=3000.0)) + + next_state = battery.step_current(state=state, current_a=12.0, dt_s=60.0) + print() + print(f"After 60 s at 12 A: SoC={next_state.soc:.3f}, DOD={next_state.dod:.3f}") + + +if __name__ == "__main__": + main() diff --git a/examples/select_motor.py b/examples/select_motor.py deleted file mode 100644 index 1f80efe3..00000000 --- a/examples/select_motor.py +++ /dev/null @@ -1,131 +0,0 @@ -"""Select commercially available brushless motors from the database. - -This script runs the powertrain co-design optimization to find the optimal theoretical -motor parameters (Kv, resistance, current limit) for a hover thrust of 4.903 N. -It then loads the brushless motor database and finds the closest matching -real-world motors sorted by internal resistance and weight. -""" - -import json -import sys -from pathlib import Path -import openmdao.api as om - -# Append project root to sys.path -sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) - -from pythrust.propellers import PropellerDatabase -from pythrust.openmdao import PropulsionComponent - - -def main(): - # 1. Load Propeller Database - db = PropellerDatabase() - db.load(Path("data/propellers/apc_202602"), strict=False) - prop_entry = db.get("APC_13x6.5E") - if prop_entry is None: - print("Error: Propeller 'APC_13x6.5E' not found.") - sys.exit(1) - - # 2. Run Co-Design Optimization to find the ideal theoretical motor - print("=== Running Theoretical Optimization ===") - prob = om.Problem() - model = prob.model - model.add_subsystem('prop', PropulsionComponent(prop_entry=prop_entry), promotes=['*']) - - model.add_design_var('kv_rpm_per_v', lower=500.0, upper=1500.0, ref=1000.0) - model.add_design_var('diameter_m', lower=0.254, upper=0.381, ref=0.3) - model.add_design_var('throttle', lower=0.2, upper=0.9, ref=0.5) - - model.add_constraint('thrust_n', equals=4.903, ref=5.0) - model.add_constraint('motor_current_a', upper=25.0, ref=25.0) - model.add_objective('battery_power_w', ref=100.0) - - prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP') - prob.driver.options['disp'] = False - prob.setup() - - # Set constants - prob.set_val('kv_rpm_per_v', 860.0) - prob.set_val('resistance_ohm', 0.0258) - prob.set_val('no_load_current_a', 1.3) - prob.set_val('current_max_a', 65.0) - prob.set_val('voltage_v', 14.8) - prob.set_val('resistance_system_ohm', 0.095) - prob.set_val('diameter_m', 0.3302) - prob.set_val('throttle', 0.5) - prob.set_val('rho', 1.225) - prob.set_val('airspeed_mps', 0.0) - - prob.run_driver() - - opt_kv = float(prob.get_val('kv_rpm_per_v')[0]) - opt_dia = float(prob.get_val('diameter_m')[0]) * 39.3701 - opt_power = float(prob.get_val('battery_power_w')[0]) - opt_throttle = float(prob.get_val('throttle')[0]) * 100.0 - opt_current = float(prob.get_val('motor_current_a')[0]) - - print(f"\nIdeal Theoretical Design:") - print(f" Target Kv : {opt_kv:.1f} RPM/V") - print(f" Target Diameter : {opt_dia:.2f} inches") - print(f" Target Hover Current: {opt_current:.2f} A") - print(f" Min Hover Power : {opt_power:.2f} W") - - # 3. Load Motor Database from PyThrust - from pythrust.motors import MotorDatabase - - motors_db_path = Path(__file__).resolve().parent.parent / "data" / "motors" - db_motors = MotorDatabase() - if not db_motors.load(motors_db_path): - print(f"\nError: Motor database not found at {motors_db_path}") - sys.exit(1) - - print(f"\n=== Loading Motor Database ({motors_db_path.name}) ===") - print(f"Loaded {db_motors.motor_count} unique reference motors.") - - # 4. Search and Filter Real Motors - min_kv = opt_kv * 0.85 - max_kv = opt_kv * 1.15 - min_max_current = opt_current * 1.5 - - results = db_motors.search( - min_kv=min_kv, - max_kv=max_kv, - min_max_current=min_max_current, - min_weight=30.0, - max_weight=180.0 - ) - - candidates = [] - for entry in results: - kv_err = abs(entry.kv - opt_kv) / opt_kv - candidates.append({ - "name": entry.name, - "manufacturer": entry.manufacturer, - "kv": entry.kv, - "resistance": entry.resistance, - "max_current": entry.max_current, - "weight_g": entry.weight_g, - "kv_error_pct": kv_err * 100.0 - }) - - # Sort candidates by winding resistance - candidates.sort(key=lambda x: (x["resistance"], x["weight_g"])) - - # 5. Display Top Matches - print(f"\n=== Top 5 Real-World Motor Matches (Ideal Kv ~ {opt_kv:.0f}) ===") - if not candidates: - print("No matching real motors found with specified criteria.") - else: - for idx, c in enumerate(candidates[:5]): - print(f"\n{idx + 1}. {c['manufacturer']} {c['name']}") - print(f" Kv : {c['kv']:.0f} RPM/V (Diff: {c['kv_error_pct']:.1f}%)") - # Convert resistance to mOhm for readability - res_mohm = c['resistance'] * 1000.0 if c['resistance'] else float('inf') - print(f" Resistance : {res_mohm:.1f} mOhm") - print(f" Weight : {c['weight_g']:.1f} g") - print(f" Max Current : {c['max_current']:.1f} A") - - -if __name__ == '__main__': - main() diff --git a/examples/select_motor_from_database.py b/examples/select_motor_from_database.py new file mode 100644 index 00000000..e9dfdd47 --- /dev/null +++ b/examples/select_motor_from_database.py @@ -0,0 +1,129 @@ +"""Select commercially available brushless motors from the database. + +This example first solves for an efficient theoretical hover design, then uses +the motor database to find real motors near the optimized Kv and current. + +Usage:: + + PYTHONPATH=. python examples/select_motor_from_database.py +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +import openmdao.api as om + +from pythrust.motors import MotorDatabase +from pythrust.openmdao import PropulsionComponent +from pythrust.propellers import PropellerDatabase + + +def main(): + # 1. Load propeller aerodynamic data. + db = PropellerDatabase() + db.load(Path("data/propellers/apc_202602"), strict=False) + prop_entry = db.get("APC_13x6.5E") + if prop_entry is None: + print("Error: Propeller 'APC_13x6.5E' not found.") + sys.exit(1) + + # 2. Run co-design optimization to find the ideal theoretical motor. + print("=== Running Theoretical Optimization ===") + prob = om.Problem() + model = prob.model + model.add_subsystem("prop", PropulsionComponent(prop_entry=prop_entry), promotes=["*"]) + + model.add_design_var("kv_rpm_per_v", lower=500.0, upper=1500.0, ref=1000.0) + model.add_design_var("diameter_m", lower=0.254, upper=0.381, ref=0.3) + model.add_design_var("throttle", lower=0.2, upper=0.9, ref=0.5) + + model.add_constraint("thrust_n", equals=4.903, ref=5.0) + model.add_constraint("motor_current_a", upper=25.0, ref=25.0) + model.add_objective("battery_power_w", ref=100.0) + + prob.driver = om.ScipyOptimizeDriver(optimizer="SLSQP") + prob.driver.options["disp"] = False + prob.setup() + + prob.set_val("kv_rpm_per_v", 860.0) + prob.set_val("resistance_ohm", 0.0258) + prob.set_val("no_load_current_a", 1.3) + prob.set_val("current_max_a", 65.0) + prob.set_val("voltage_v", 14.8) + prob.set_val("resistance_system_ohm", 0.095) + prob.set_val("diameter_m", 0.3302) + prob.set_val("throttle", 0.5) + prob.set_val("rho", 1.225) + prob.set_val("airspeed_mps", 0.0) + + prob.run_driver() + + opt_kv = float(prob.get_val("kv_rpm_per_v")[0]) + opt_dia = float(prob.get_val("diameter_m")[0]) * 39.3701 + opt_power = float(prob.get_val("battery_power_w")[0]) + opt_current = float(prob.get_val("motor_current_a")[0]) + + print("\nIdeal Theoretical Design:") + print(f" Target Kv : {opt_kv:.1f} RPM/V") + print(f" Target Diameter : {opt_dia:.2f} inches") + print(f" Target Hover Current: {opt_current:.2f} A") + print(f" Min Hover Power : {opt_power:.2f} W") + + # 3. Load motor database and find catalog candidates near the optimum. + motors_db_path = Path(__file__).resolve().parent.parent / "data" / "motors" + db_motors = MotorDatabase() + if not db_motors.load(motors_db_path): + print(f"\nError: Motor database not found at {motors_db_path}") + sys.exit(1) + + print(f"\n=== Loading Motor Database ({motors_db_path.name}) ===") + print(f"Loaded {db_motors.motor_count} unique reference motors.") + + min_kv = opt_kv * 0.85 + max_kv = opt_kv * 1.15 + min_max_current = opt_current * 1.5 + + results = db_motors.search( + min_kv=min_kv, + max_kv=max_kv, + min_max_current=min_max_current, + min_weight=30.0, + max_weight=180.0, + ) + + candidates = [] + for entry in results: + kv_err = abs(entry.kv - opt_kv) / opt_kv + candidates.append( + { + "name": entry.name, + "manufacturer": entry.manufacturer, + "kv": entry.kv, + "resistance": entry.resistance, + "max_current": entry.max_current, + "weight_g": entry.weight_g, + "kv_error_pct": kv_err * 100.0, + } + ) + + candidates.sort(key=lambda x: (x["resistance"], x["weight_g"])) + + print(f"\n=== Top 5 Real-World Motor Matches (Ideal Kv ~ {opt_kv:.0f}) ===") + if not candidates: + print("No matching real motors found with specified criteria.") + else: + for idx, c in enumerate(candidates[:5]): + print(f"\n{idx + 1}. {c['manufacturer']} {c['name']}") + print(f" Kv : {c['kv']:.0f} RPM/V (Diff: {c['kv_error_pct']:.1f}%)") + res_mohm = c["resistance"] * 1000.0 if c["resistance"] else float("inf") + print(f" Resistance : {res_mohm:.1f} mOhm") + print(f" Weight : {c['weight_g']:.1f} g") + print(f" Max Current : {c['max_current']:.1f} A") + + +if __name__ == "__main__": + main() diff --git a/examples/simulate_pybamm_mission.py b/examples/simulate_pybamm_mission.py deleted file mode 100644 index 62df5b63..00000000 --- a/examples/simulate_pybamm_mission.py +++ /dev/null @@ -1,197 +0,0 @@ -import matplotlib -matplotlib.use('Agg') - -import matplotlib.pyplot as plt -import numpy as np -from pathlib import Path -import sys -import pybamm - -# Add project root to sys.path -sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) - -from pythrust.propellers import PropellerDatabase -from pythrust.propulsion import ( - BatterySpec, - MotorSpec, - PropellerSpec, - SystemSpec, - PropulsionSolver, -) - -def get_mission_throttle(t): - """Define a simplified UAV flight mission throttle profile. - - - 0 to 60s (1 min) : Takeoff & Climb (75% throttle) - - 60s to 660s (10 min) : Cruise (50% throttle) - - 660s to 720s (1 min) : Descend & Landing (35% throttle) - """ - if t < 60.0: - return 0.75 - elif t < 660.0: - return 0.50 - elif t < 720.0: - return 0.35 - else: - return 0.0 - -def main(): - # 1. Load PyThrust databases - db = PropellerDatabase() - db.load(Path("data/propellers/apc_202602")) - prop_entry = db.get("APC_13x6.5E") - - motor = MotorSpec(kv_rpm_per_v=860.0, resistance_ohm=0.0258, no_load_current_a=1.3, current_max_a=65.0) - system = SystemSpec(resistance_ohm=0.05) - propeller = PropellerSpec(diameter_m=0.3302) - solver = PropulsionSolver() - - # 2. Setup PyBaMM Lithium-ion Model (Single Particle Model) - print("Initializing PyBaMM Single Particle Model (SPM)...") - model = pybamm.lithium_ion.SPM() - parameter_values = model.default_parameter_values - - # Set the current as an input parameter for dynamic stepping - parameter_values.update({"Current function [A]": "[input]"}) - - sim = pybamm.Simulation(model, parameter_values=parameter_values) - - # Battery configuration scaling details - cells_series = 4 - capacity_reference = parameter_values["Nominal cell capacity [A.h]"] # ~0.68 Ah - capacity_pack = 5.0 # 5000 mAh target pack - current_scaling_ratio = capacity_reference / capacity_pack - - # Simulation settings - dt = 1.0 # 1s step - total_time = 720.0 - time_steps = np.arange(0.0, total_time, dt) - - times = [] - throttles = [] - voltages = [] - currents = [] - thrusts_g = [] - soc_history = [] - - # Initial states - battery_voltage = 4.2 * cells_series # starting voltage - last_current = 0.0 - - print("Running dynamic mission profile simulation using PyBaMM...") - - for t in time_steps: - throttle = get_mission_throttle(t) - if throttle <= 0.0: - break - - battery_spec = BatterySpec(voltage_v=battery_voltage) - - # Solve operating point at this voltage and throttle - pt = solver.solve_operating_point( - motor=motor, - battery=battery_spec, - system=system, - propeller=propeller, - prop_entry=prop_entry, - rho=1.225, - airspeed_mps=0.0, - throttle=throttle - ) - - if not pt.is_feasible: - print(f"System infeasible at t = {t}s due to: {pt.infeasible_reason}") - break - - current = pt.battery_power_w / battery_voltage - - # Scale current for the reference cell in PyBaMM - current_cell = current * current_scaling_ratio - - # Step PyBaMM model - try: - sim.step(dt=dt, inputs={"Current function [A]": current_cell}) - except Exception as e: - print(f"PyBaMM step failed (battery depleted or voltage limit reached) at t = {t}s: {e}") - break - - # Retrieve state from solver solution - v_cell = sim.solution["Terminal voltage [V]"].data[-1] - battery_voltage = v_cell * cells_series - - discharged_ah = sim.solution["Discharge capacity [A.h]"].data[-1] - soc = 1.0 - (discharged_ah / capacity_reference) - - times.append(t) - throttles.append(throttle * 100.0) - voltages.append(battery_voltage) - currents.append(current) - thrusts_g.append(pt.thrust_n * 1000.0 / 9.80665) - soc_history.append(soc * 100.0) - - # Stop if SoC drops below 2% - if soc <= 0.02: - print(f"Battery depleted (2% SoC) at t = {t}s") - break - - # Convert lists to arrays - times = np.array(times) / 60.0 # to minutes - throttles = np.array(throttles) - voltages = np.array(voltages) - currents = np.array(currents) - thrusts_g = np.array(thrusts_g) - soc_history = np.array(soc_history) - - # Plotting - fig, axs = plt.subplots(2, 2, figsize=(12, 10)) - - # Subplot 1 (Top-Left): Throttle Profile & Voltage - color = 'C0' - axs[0, 0].plot(times, throttles, color=color, label='Throttle (%)', linewidth=2) - axs[0, 0].set_ylabel('Throttle (%)', color=color) - axs[0, 0].tick_params(axis='y', labelcolor=color) - axs[0, 0].grid(True) - axs[0, 0].set_xlabel('Time (minutes)') - axs[0, 0].set_title('Throttle & Voltage') - - ax0_twin = axs[0, 0].twinx() - color = 'C3' - ax0_twin.plot(times, voltages, color=color, label='Voltage (V)', linestyle='--', linewidth=1.8) - ax0_twin.set_ylabel('Terminal Voltage (V)', color=color) - ax0_twin.tick_params(axis='y', labelcolor=color) - - # Subplot 2 (Top-Right): State of Charge (%) - axs[0, 1].plot(times, soc_history, color='C4', linewidth=2) - axs[0, 1].set_ylabel('State of Charge (%)') - axs[0, 1].set_xlabel('Time (minutes)') - axs[0, 1].grid(True) - axs[0, 1].set_title('State of Charge (%)') - - # Subplot 3 (Bottom-Left): Battery Current Draw (A) - axs[1, 0].plot(times, currents, color='C1', linewidth=2) - axs[1, 0].set_ylabel('Battery Current (A)') - axs[1, 0].set_xlabel('Time (minutes)') - axs[1, 0].grid(True) - axs[1, 0].set_title('Battery Current Draw') - - # Subplot 4 (Bottom-Right): Produced Thrust (g) - axs[1, 1].plot(times, thrusts_g, color='C2', linewidth=2) - axs[1, 1].set_ylabel('Thrust (g)') - axs[1, 1].set_xlabel('Time (minutes)') - axs[1, 1].grid(True) - axs[1, 1].set_title('Produced Thrust') - - fig.suptitle('UAV Propulsion System Simulation: PyBaMM Electrochemical Battery Model', fontsize=14, y=0.98) - plt.tight_layout() - - # Save the output image - output_dir = Path("docs/images") - output_dir.mkdir(parents=True, exist_ok=True) - output_image = output_dir / "pybamm_mission_results.png" - plt.savefig(output_image, dpi=150, bbox_inches='tight') - plt.close() - - print(f"Successfully generated and saved PyBaMM simulation plot to: {output_image.resolve()}") - -if __name__ == '__main__': - main() diff --git a/mkdocs.yml b/mkdocs.yml index da791798..7f3c963b 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -39,6 +39,8 @@ extra_css: - stylesheets/extra.css markdown_extensions: + - attr_list + - md_in_html - admonition - pymdownx.details - pymdownx.superfences: @@ -70,7 +72,8 @@ nav: - Examples: examples.md - Reference: - API Reference: api_reference.md - - Mathematical Model: theory.md + - Propulsion and Battery Theory: theory.md + - Battery Model: battery_model.md - Component Databases: databases.md - Developer Guide: - Development & Testing: development.md diff --git a/pythrust/battery/__init__.py b/pythrust/battery/__init__.py new file mode 100644 index 00000000..663d1e96 --- /dev/null +++ b/pythrust/battery/__init__.py @@ -0,0 +1,12 @@ +"""Battery models for propulsion and mission analysis.""" + +from .fixed import FixedVoltageBattery +from .rate_map import RateMapBattery +from .state import BatteryPoint, BatteryState + +__all__ = [ + "BatteryPoint", + "BatteryState", + "FixedVoltageBattery", + "RateMapBattery", +] diff --git a/pythrust/battery/fixed.py b/pythrust/battery/fixed.py new file mode 100644 index 00000000..a8063726 --- /dev/null +++ b/pythrust/battery/fixed.py @@ -0,0 +1,34 @@ +"""Fixed-voltage battery model.""" + +from __future__ import annotations + +from dataclasses import dataclass + +from .state import BatteryState + + +@dataclass(frozen=True) +class FixedVoltageBattery: + """Battery model with constant pack voltage. + + This is the historical PyThrust battery model under a more explicit name. + """ + + voltage_v: float + discharge_efficiency: float = 1.0 + + def terminal_voltage( + self, + current_a: float = 0.0, + state: BatteryState | None = None, + ) -> float: + """Return pack voltage for the requested current and state.""" + return self.voltage_v + + def terminal_power( + self, + current_a: float, + state: BatteryState | None = None, + ) -> float: + """Return battery-side power draw using the fixed-voltage efficiency.""" + return self.voltage_v * current_a / max(1e-12, self.discharge_efficiency) diff --git a/pythrust/battery/rate_map.py b/pythrust/battery/rate_map.py new file mode 100644 index 00000000..09c28674 --- /dev/null +++ b/pythrust/battery/rate_map.py @@ -0,0 +1,255 @@ +"""Rate-map battery model.""" + +from __future__ import annotations + +from dataclasses import dataclass +import json +import math +from pathlib import Path +from typing import Sequence + +import numpy as np + +from .state import BatteryPoint, BatteryState + + +@dataclass(frozen=True) +class RateMapBattery: + """Equivalent-circuit battery model driven by OCV and resistance curves.""" + + name: str + capacity_ah: float + cutoff_voltage_v: float + charge_voltage_v: float + max_current_a: float + series: int + parallel: int + dod: Sequence[float] + ocv_v: Sequence[float] + resistance_ohm: Sequence[float] + source: str | None = None + + def __post_init__(self) -> None: + dod = tuple(float(v) for v in self.dod) + ocv_v = tuple(float(v) for v in self.ocv_v) + resistance_ohm = tuple(float(v) for v in self.resistance_ohm) + + if len(dod) < 2: + raise ValueError("dod must contain at least two points") + if not (len(dod) == len(ocv_v) == len(resistance_ohm)): + raise ValueError("dod, ocv_v, and resistance_ohm must have the same length") + if any(v < 0.0 or v > 1.0 for v in dod): + raise ValueError("dod values must be between 0 and 1") + if any(b <= a for a, b in zip(dod, dod[1:])): + raise ValueError("dod values must be strictly increasing") + if any(v <= 0.0 for v in resistance_ohm): + raise ValueError("resistance_ohm values must be positive") + if self.capacity_ah <= 0.0: + raise ValueError("capacity_ah must be positive") + if self.max_current_a <= 0.0: + raise ValueError("max_current_a must be positive") + if self.series <= 0 or self.parallel <= 0: + raise ValueError("series and parallel counts must be positive") + + object.__setattr__(self, "dod", dod) + object.__setattr__(self, "ocv_v", ocv_v) + object.__setattr__(self, "resistance_ohm", resistance_ohm) + + @classmethod + def from_json( + cls, + path: str | Path, + *, + series: int = 1, + parallel: int = 1, + ) -> "RateMapBattery": + """Load a cell dataset and apply pack series/parallel topology.""" + with Path(path).open("r", encoding="utf-8") as f: + data = json.load(f) + + cell = data["cell"] + curves = data["curves"] + return cls( + name=data.get("name", "Rate-map cell"), + source=data.get("source"), + capacity_ah=cell["capacity_ah"], + cutoff_voltage_v=cell["cutoff_voltage_v"], + charge_voltage_v=cell["charge_voltage_v"], + max_current_a=cell["max_current_a"], + series=series, + parallel=parallel, + dod=curves["dod"], + ocv_v=curves["ocv_v"], + resistance_ohm=curves["resistance_ohm"], + ) + + @property + def capacity_as(self) -> float: + """Cell capacity in ampere-seconds.""" + return self.capacity_ah * 3600.0 + + @property + def rated_current_a(self) -> float: + """Cell 1C current in amps.""" + return self.capacity_ah + + def ocv(self, dod: float) -> float: + """Interpolate open-circuit cell voltage at depth of discharge.""" + return float(np.interp(self._clip_dod(dod), self.dod, self.ocv_v)) + + def resistance(self, dod: float) -> float: + """Interpolate cell resistance at depth of discharge.""" + return float(np.interp(self._clip_dod(dod), self.dod, self.resistance_ohm)) + + def terminal_voltage(self, current_a: float, state: BatteryState) -> float: + """Return pack terminal voltage at pack current and battery state.""" + return self.state_at_current(state=state, current_a=current_a).terminal_voltage_v + + def terminal_power(self, current_a: float, state: BatteryState) -> float: + """Return pack terminal power at pack current and battery state.""" + return self.state_at_current(state=state, current_a=current_a).power_w + + def state_at_current(self, state: BatteryState, current_a: float) -> BatteryPoint: + """Evaluate battery point state for a specified pack current.""" + cell_current_a = current_a / self.parallel + dod = state.dod + ocv_v = self.ocv(dod) + resistance_ohm = self.resistance(dod) + cell_voltage_v = ocv_v - resistance_ohm * cell_current_a + return self._point_from_cell_state( + state=state, + cell_current_a=cell_current_a, + cell_voltage_v=cell_voltage_v, + ocv_v=ocv_v, + resistance_ohm=resistance_ohm, + ) + + def state_at_c_rate(self, state: BatteryState, c_rate: float) -> BatteryPoint: + """Evaluate battery point state for a specified cell C-rate.""" + return self.state_at_current( + state=state, + current_a=c_rate * self.rated_current_a * self.parallel, + ) + + def state_at_voltage(self, state: BatteryState, voltage_v: float) -> BatteryPoint: + """Evaluate battery point state for a specified pack voltage.""" + cell_voltage_v = voltage_v / self.series + dod = state.dod + ocv_v = self.ocv(dod) + resistance_ohm = self.resistance(dod) + cell_current_a = (ocv_v - cell_voltage_v) / resistance_ohm + return self._point_from_cell_state( + state=state, + cell_current_a=cell_current_a, + cell_voltage_v=cell_voltage_v, + ocv_v=ocv_v, + resistance_ohm=resistance_ohm, + ) + + def state_at_power(self, state: BatteryState, power_w: float) -> BatteryPoint: + """Evaluate battery point state for a specified pack power draw.""" + dod = state.dod + ocv_v = self.ocv(dod) + resistance_ohm = self.resistance(dod) + cell_power_w = power_w / (self.series * self.parallel) + discriminant = ocv_v**2 - 4.0 * resistance_ohm * cell_power_w + + if discriminant < 0.0: + return BatteryPoint( + terminal_voltage_v=0.0, + current_a=0.0, + power_w=0.0, + cell_voltage_v=0.0, + cell_current_a=0.0, + c_rate=0.0, + efficiency=0.0, + ocv_v=ocv_v, + resistance_ohm=resistance_ohm, + is_feasible=False, + infeasible_reason="power_limit", + ) + + cell_current_a = (ocv_v - math.sqrt(discriminant)) / (2.0 * resistance_ohm) + cell_voltage_v = ocv_v - resistance_ohm * cell_current_a + return self._point_from_cell_state( + state=state, + cell_current_a=cell_current_a, + cell_voltage_v=cell_voltage_v, + ocv_v=ocv_v, + resistance_ohm=resistance_ohm, + ) + + def state_at_load_resistance(self, state: BatteryState, resistance_ohm: float) -> BatteryPoint: + """Evaluate battery point state for a specified pack load resistance.""" + if resistance_ohm <= 0.0: + raise ValueError("load resistance must be positive") + + dod = state.dod + ocv_v = self.ocv(dod) + cell_internal_ohm = self.resistance(dod) + cell_load_ohm = resistance_ohm * self.parallel / self.series + cell_current_a = ocv_v / (cell_internal_ohm + cell_load_ohm) + cell_voltage_v = ocv_v - cell_internal_ohm * cell_current_a + return self._point_from_cell_state( + state=state, + cell_current_a=cell_current_a, + cell_voltage_v=cell_voltage_v, + ocv_v=ocv_v, + resistance_ohm=cell_internal_ohm, + ) + + def step_current(self, state: BatteryState, current_a: float, dt_s: float) -> BatteryState: + """Advance state using constant pack current over a time step.""" + if dt_s < 0.0: + raise ValueError("dt_s must be non-negative") + cell_current_a = current_a / self.parallel + dod_next = state.dod + cell_current_a * dt_s / self.capacity_as + return BatteryState.from_dod(min(1.0, max(0.0, dod_next))) + + def step_power(self, state: BatteryState, power_w: float, dt_s: float) -> BatteryState: + """Advance state using constant pack power over a time step.""" + point = self.state_at_power(state=state, power_w=power_w) + if not point.is_feasible: + raise ValueError(point.infeasible_reason or "infeasible battery power") + return self.step_current(state=state, current_a=point.current_a, dt_s=dt_s) + + def _point_from_cell_state( + self, + state: BatteryState, + cell_current_a: float, + cell_voltage_v: float, + ocv_v: float, + resistance_ohm: float, + ) -> BatteryPoint: + pack_current_a = cell_current_a * self.parallel + pack_voltage_v = cell_voltage_v * self.series + c_rate = cell_current_a / self.rated_current_a + efficiency = cell_voltage_v / ocv_v if ocv_v > 0.0 else 0.0 + reason = None + + if cell_current_a > self.max_current_a: + reason = "current_limit" + if reason is None and cell_voltage_v < self.cutoff_voltage_v: + reason = "voltage_cutoff" + if reason is None and cell_voltage_v > self.charge_voltage_v: + reason = "voltage_limit" + if reason is None and (state.dod < 0.0 or state.dod > 1.0): + reason = "state_limit" + + return BatteryPoint( + terminal_voltage_v=float(pack_voltage_v), + current_a=float(pack_current_a), + power_w=float(pack_voltage_v * pack_current_a), + cell_voltage_v=float(cell_voltage_v), + cell_current_a=float(cell_current_a), + c_rate=float(c_rate), + efficiency=float(efficiency), + ocv_v=float(ocv_v), + resistance_ohm=float(resistance_ohm), + is_feasible=reason is None, + infeasible_reason=reason, + ) + + @staticmethod + def _clip_dod(dod: float) -> float: + return min(1.0, max(0.0, float(dod))) diff --git a/pythrust/battery/state.py b/pythrust/battery/state.py new file mode 100644 index 00000000..778b534e --- /dev/null +++ b/pythrust/battery/state.py @@ -0,0 +1,47 @@ +"""Battery state and point-result data structures.""" + +from __future__ import annotations + +from dataclasses import dataclass +import math +from typing import Optional + + +@dataclass(frozen=True) +class BatteryState: + """Battery state for mission or point-performance analysis.""" + + soc: float = 1.0 + + def __post_init__(self) -> None: + if not math.isfinite(self.soc) or self.soc < 0.0 or self.soc > 1.0: + raise ValueError("soc must be between 0 and 1") + + @classmethod + def from_dod(cls, dod: float) -> "BatteryState": + """Build a battery state from depth of discharge.""" + if not math.isfinite(dod) or dod < 0.0 or dod > 1.0: + raise ValueError("dod must be between 0 and 1") + return cls(soc=1.0 - dod) + + @property + def dod(self) -> float: + """Depth of discharge, from 0 full to 1 empty.""" + return 1.0 - self.soc + + +@dataclass(frozen=True) +class BatteryPoint: + """Battery point-state result.""" + + terminal_voltage_v: float + current_a: float + power_w: float + cell_voltage_v: float + cell_current_a: float + c_rate: float + efficiency: float + ocv_v: float + resistance_ohm: float + is_feasible: bool + infeasible_reason: Optional[str] = None diff --git a/pythrust/openmdao/components.py b/pythrust/openmdao/components.py index cbc6534c..21609404 100644 --- a/pythrust/openmdao/components.py +++ b/pythrust/openmdao/components.py @@ -96,7 +96,7 @@ def compute(self, inputs: om.Vector, outputs: om.Vector) -> None: outputs['rpm'] = op.rpm outputs['thrust_n'] = op.thrust_n outputs['torque_nm'] = op.torque_nm - outputs['battery_current_a'] = op.battery_power_w / max(1e-6, battery.voltage_v) + outputs['battery_current_a'] = op.battery_current_a outputs['battery_power_w'] = op.battery_power_w outputs['motor_current_a'] = op.motor_current_a outputs['motor_voltage_v'] = op.motor_voltage_v diff --git a/pythrust/propulsion/models.py b/pythrust/propulsion/models.py deleted file mode 100644 index 11b9270a..00000000 --- a/pythrust/propulsion/models.py +++ /dev/null @@ -1,124 +0,0 @@ -"""Core propulsion model data structures.""" - -from __future__ import annotations - -from dataclasses import dataclass -from typing import Optional - - -import math - - -@dataclass(frozen=True) -class MotorSpec: - """Motor electrical parameters (supporting both 1st and 2nd order models). - - All values are taken directly from the manufacturer datasheet or calibrated. - - Units: - - kv_rpm_per_v: RPM / V - - resistance_ohm: ohm - - no_load_current_a: A - - current_max_a: A - """ - - kv_rpm_per_v: float - resistance_ohm: float - no_load_current_a: float - current_max_a: float - - # Optional 2nd-order parameters (Drela / QPROP model) - torque_constant_kv_ratio: float = 1.0 # Kq / Kv ratio (default: 1.0) - magnetic_lag_tau: float = 0.0 # tau (magnetic lag time constant) - no_load_current_linear: float = 0.0 # Io1 (A/(rad/s)) - no_load_current_quadratic: float = 0.0 # Io2 (A/(rad/s)^2) - resistance_quadratic: float = 0.0 # R2 (Ohms/Amp^2) - - # Simplified power-law iron loss parameters - no_load_voltage_v: float = 10.0 # Datasheet measurement voltage - iron_loss_exponent: float = 0.0 # If >0, scales I0 by (RPM / RPM_0)^exponent - - def get_no_load_current(self, rpm: float) -> float: - """Calculate the no-load current at a given shaft speed (RPM).""" - if rpm <= 0.0: - return self.no_load_current_a - - # 1. Check if power-law exponent is specified - if self.iron_loss_exponent > 0.0: - rpm_0 = self.kv_rpm_per_v * self.no_load_voltage_v - if rpm_0 > 0.0: - return self.no_load_current_a * (rpm / rpm_0) ** self.iron_loss_exponent - - # 2. Otherwise use Drela's quadratic speed model: Io0 + Io1*omega + Io2*omega^2 - omega = rpm * (math.pi / 30.0) - return ( - self.no_load_current_a - + self.no_load_current_linear * omega - + self.no_load_current_quadratic * (omega ** 2) - ) - - def get_winding_resistance(self, current_a: float) -> float: - """Calculate winding resistance including current-dependent heating.""" - if self.resistance_quadratic <= 0.0: - return self.resistance_ohm - return self.resistance_ohm + self.resistance_quadratic * (current_a ** 2) - - -@dataclass(frozen=True) -class BatterySpec: - """Battery pack parameters. - - Units: - - voltage_v: V - - discharge_efficiency: 0-1 - """ - - voltage_v: float - discharge_efficiency: float = 1.0 - - -@dataclass(frozen=True) -class SystemSpec: - """System transmission/line electrical parameters. - - Units: - - resistance_ohm: ohm - """ - - resistance_ohm: float = 0.0 - - -@dataclass(frozen=True) -class PropellerSpec: - """Propeller geometry. - - Units: - - diameter_m: m - - pitch_m: m (optional) - """ - - diameter_m: float - blade_count: int = 2 - pitch_m: Optional[float] = None - - -@dataclass(frozen=True) -class OperatingPoint: - """Solved operating point for a given condition.""" - - rpm: float - advance_ratio: float - ct: float - cp: float - thrust_n: float - torque_nm: float - shaft_power_w: float - motor_power_w: float - battery_power_w: float - motor_current_a: float - 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/models/__init__.py b/pythrust/propulsion/models/__init__.py new file mode 100644 index 00000000..7af2ee7e --- /dev/null +++ b/pythrust/propulsion/models/__init__.py @@ -0,0 +1,15 @@ +"""Propulsion model data structures.""" + +from .battery import BatterySpec +from .motor import MotorSpec +from .operating_point import OperatingPoint +from .propeller import PropellerSpec +from .system import SystemSpec + +__all__ = [ + "BatterySpec", + "MotorSpec", + "OperatingPoint", + "PropellerSpec", + "SystemSpec", +] diff --git a/pythrust/propulsion/models/battery.py b/pythrust/propulsion/models/battery.py new file mode 100644 index 00000000..cdbc8296 --- /dev/null +++ b/pythrust/propulsion/models/battery.py @@ -0,0 +1,7 @@ +"""Compatibility battery model exports for propulsion.""" + +from pythrust.battery import FixedVoltageBattery + +BatterySpec = FixedVoltageBattery + +__all__ = ["BatterySpec"] diff --git a/pythrust/propulsion/models/motor.py b/pythrust/propulsion/models/motor.py new file mode 100644 index 00000000..0ccf1a11 --- /dev/null +++ b/pythrust/propulsion/models/motor.py @@ -0,0 +1,59 @@ +"""Motor model data structures.""" + +from __future__ import annotations + +from dataclasses import dataclass +import math + + +@dataclass(frozen=True) +class MotorSpec: + """Motor electrical parameters (supporting both 1st and 2nd order models). + + All values are taken directly from the manufacturer datasheet or calibrated. + + Units: + - kv_rpm_per_v: RPM / V + - resistance_ohm: ohm + - no_load_current_a: A + - current_max_a: A + """ + + kv_rpm_per_v: float + resistance_ohm: float + no_load_current_a: float + current_max_a: float + + # Optional 2nd-order parameters (Drela / QPROP model) + torque_constant_kv_ratio: float = 1.0 + magnetic_lag_tau: float = 0.0 + no_load_current_linear: float = 0.0 + no_load_current_quadratic: float = 0.0 + resistance_quadratic: float = 0.0 + + # Simplified power-law iron loss parameters + no_load_voltage_v: float = 10.0 + iron_loss_exponent: float = 0.0 + + def get_no_load_current(self, rpm: float) -> float: + """Calculate the no-load current at a given shaft speed (RPM).""" + if rpm <= 0.0: + return self.no_load_current_a + + if self.iron_loss_exponent > 0.0: + rpm_0 = self.kv_rpm_per_v * self.no_load_voltage_v + if rpm_0 > 0.0: + return self.no_load_current_a * (rpm / rpm_0) ** self.iron_loss_exponent + + omega = rpm * (math.pi / 30.0) + return ( + self.no_load_current_a + + self.no_load_current_linear * omega + + self.no_load_current_quadratic * (omega**2) + ) + + def get_winding_resistance(self, current_a: float) -> float: + """Calculate winding resistance including current-dependent heating.""" + if self.resistance_quadratic <= 0.0: + return self.resistance_ohm + return self.resistance_ohm + self.resistance_quadratic * (current_a**2) diff --git a/pythrust/propulsion/models/operating_point.py b/pythrust/propulsion/models/operating_point.py new file mode 100644 index 00000000..a1ad225d --- /dev/null +++ b/pythrust/propulsion/models/operating_point.py @@ -0,0 +1,32 @@ +"""Solved propulsion operating point data structures.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Optional + + +@dataclass(frozen=True) +class OperatingPoint: + """Solved operating point for a given condition.""" + + rpm: float + advance_ratio: float + ct: float + cp: float + thrust_n: float + torque_nm: float + shaft_power_w: float + motor_power_w: float + battery_power_w: float + motor_current_a: float + 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 + battery_voltage_v: float = 0.0 + battery_current_a: float = 0.0 + battery_c_rate: float = 0.0 + battery_efficiency: float = 1.0 diff --git a/pythrust/propulsion/models/propeller.py b/pythrust/propulsion/models/propeller.py new file mode 100644 index 00000000..f15869b7 --- /dev/null +++ b/pythrust/propulsion/models/propeller.py @@ -0,0 +1,20 @@ +"""Propeller model data structures.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Optional + + +@dataclass(frozen=True) +class PropellerSpec: + """Propeller geometry. + + Units: + - diameter_m: m + - pitch_m: m (optional) + """ + + diameter_m: float + blade_count: int = 2 + pitch_m: Optional[float] = None diff --git a/pythrust/propulsion/models/system.py b/pythrust/propulsion/models/system.py new file mode 100644 index 00000000..3e2928e2 --- /dev/null +++ b/pythrust/propulsion/models/system.py @@ -0,0 +1,16 @@ +"""Electrical system model data structures.""" + +from __future__ import annotations + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class SystemSpec: + """System transmission/line electrical parameters. + + Units: + - resistance_ohm: ohm + """ + + resistance_ohm: float = 0.0 diff --git a/pythrust/propulsion/solver.py b/pythrust/propulsion/solver.py index 9b844c03..c8be3694 100644 --- a/pythrust/propulsion/solver.py +++ b/pythrust/propulsion/solver.py @@ -9,6 +9,7 @@ from scipy.optimize import root_scalar +from pythrust.battery import BatteryPoint, BatteryState from pythrust.propellers.database import PropellerEntry from .models import BatterySpec, MotorSpec, OperatingPoint, PropellerSpec, SystemSpec @@ -63,6 +64,7 @@ def solve_operating_point( rho: float, airspeed_mps: float, throttle: float, + battery_state: BatteryState | None = None, ) -> OperatingPoint: if throttle <= 0.0: return OperatingPoint( @@ -79,8 +81,14 @@ def solve_operating_point( motor_voltage_v=0.0, is_feasible=False, infeasible_reason="throttle<=0", + battery_voltage_v=self._battery_voltage(battery, 0.0, battery_state), + battery_current_a=0.0, + battery_c_rate=0.0, + battery_efficiency=self._battery_efficiency(battery, None), ) + self._validate_battery_state(battery, battery_state) + def g(rpm: float) -> float: ct, cp, j, torque_nm, current_a, v_back = self._evaluate_state( motor, @@ -93,7 +101,9 @@ def g(rpm: float) -> float: if cp <= 0.0 or ct < 0.0 or j < 0.0: return float("inf") v_motor = v_back + current_a * motor.get_winding_resistance(current_a) - v_applied = throttle * battery.voltage_v + battery_current_for_voltage_a = throttle * current_a + battery_voltage_v = self._battery_voltage(battery, battery_current_for_voltage_a, battery_state) + v_applied = throttle * battery_voltage_v return v_motor + current_a * system.resistance_ohm - v_applied rpm_min = max(self.config.rpm_min, 1.0) @@ -101,7 +111,7 @@ def g(rpm: float) -> float: j_max = self._estimate_j_max(prop_entry) rpm_min = max(rpm_min, airspeed_mps / (j_max * propeller.diameter_m) * 60.0) - rpm_max = motor.kv_rpm_per_v * battery.voltage_v * throttle + rpm_max = motor.kv_rpm_per_v * self._battery_voltage(battery, 0.0, battery_state) * throttle rpm_max = max(rpm_max, rpm_min * 1.2) * self.config.rpm_max_margin g_min = g(rpm_min) @@ -118,6 +128,8 @@ def g(rpm: float) -> float: airspeed_mps, rpm_min if abs(g_min) < abs(g_max) else rpm_max, reason="no_bracket", + battery_state=battery_state, + throttle=throttle, ) result = root_scalar( @@ -140,6 +152,8 @@ def g(rpm: float) -> float: airspeed_mps, result.root, reason="no_convergence", + battery_state=battery_state, + throttle=throttle, ) return self._build_point( @@ -151,6 +165,8 @@ def g(rpm: float) -> float: rho, airspeed_mps, result.root, + battery_state=battery_state, + throttle=throttle, ) def _evaluate_state( @@ -174,6 +190,8 @@ def _build_point( rho: float, airspeed_mps: float, rpm: float, + battery_state: BatteryState | None = None, + throttle: float = 1.0, ) -> OperatingPoint: ct, cp, j, torque_nm, current_a, v_back = self._evaluate_state( motor, @@ -189,7 +207,27 @@ def _build_point( shaft_power_w = cp * rho * (n**3) * (propeller.diameter_m**5) motor_voltage_v = v_back + current_a * motor.get_winding_resistance(current_a) 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) + battery_power_w = self._battery_power( + battery, + current_a=current_a, + motor_power_w=motor_power_w, + system=system, + ) + battery_current_for_voltage_a = throttle * current_a + battery_point = self._battery_point(battery, battery_current_for_voltage_a, battery_state) + battery_voltage_v = self._battery_voltage_from_point( + battery, + battery_current_for_voltage_a, + battery_state, + battery_point, + ) + battery_current_a = self._battery_current_from_point_or_power( + battery_power_w, + battery_voltage_v, + battery_point, + ) + battery_c_rate = battery_point.c_rate if battery_point is not None else 0.0 + battery_efficiency = self._battery_efficiency(battery, battery_point) # Efficiency calculations if shaft_power_w > 0.0: @@ -207,18 +245,18 @@ def _build_point( else: system_efficiency = 0.0 - feasible = True reason = None if current_a > motor.current_max_a: - feasible = False reason = "current_limit" - if ct < 0.0 or cp < 0.0 or j < 0.0: - feasible = False + if reason is None and battery_point is not None and not battery_point.is_feasible: + reason = f"battery_{battery_point.infeasible_reason}" + if reason is None and (ct < 0.0 or cp < 0.0 or j < 0.0): reason = "invalid_coefficients" - if (propeller_efficiency > 1.0001 or propeller_efficiency < 0.0 or + if reason is None and ( + 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 + system_efficiency > 1.0001 or system_efficiency < 0.0 + ): reason = "invalid_efficiency" return OperatingPoint( @@ -233,11 +271,15 @@ def _build_point( battery_power_w=float(battery_power_w), motor_current_a=float(current_a), motor_voltage_v=float(motor_voltage_v), - is_feasible=feasible, + is_feasible=reason is None, infeasible_reason=reason, propeller_efficiency=float(propeller_efficiency), motor_efficiency=float(motor_efficiency), system_efficiency=float(system_efficiency), + battery_voltage_v=float(battery_voltage_v), + battery_current_a=float(battery_current_a), + battery_c_rate=float(battery_c_rate), + battery_efficiency=float(battery_efficiency), ) def _build_infeasible_point( @@ -251,6 +293,8 @@ def _build_infeasible_point( airspeed_mps: float, rpm: float, reason: str, + battery_state: BatteryState | None = None, + throttle: float = 1.0, ) -> OperatingPoint: point = self._build_point( motor, @@ -261,6 +305,8 @@ def _build_infeasible_point( rho, airspeed_mps, rpm, + battery_state=battery_state, + throttle=throttle, ) return OperatingPoint( rpm=point.rpm, @@ -279,8 +325,81 @@ def _build_infeasible_point( propeller_efficiency=point.propeller_efficiency, motor_efficiency=point.motor_efficiency, system_efficiency=point.system_efficiency, + battery_voltage_v=point.battery_voltage_v, + battery_current_a=point.battery_current_a, + battery_c_rate=point.battery_c_rate, + battery_efficiency=point.battery_efficiency, ) + @staticmethod + def _validate_battery_state(battery: BatterySpec, state: BatteryState | None) -> None: + if state is None and not hasattr(battery, "voltage_v"): + raise ValueError("battery_state is required for dynamic battery models") + + def _battery_voltage( + self, + battery: BatterySpec, + current_a: float, + state: BatteryState | None, + ) -> float: + point = self._battery_point(battery, current_a, state) + if point is not None: + return point.terminal_voltage_v + return float(battery.terminal_voltage(current_a=current_a, state=state)) + + def _battery_voltage_from_point( + self, + battery: BatterySpec, + current_a: float, + state: BatteryState | None, + point: BatteryPoint | None, + ) -> float: + if point is not None: + return point.terminal_voltage_v + return self._battery_voltage(battery, current_a, state) + + @staticmethod + def _battery_current_from_point_or_power( + battery_power_w: float, + battery_voltage_v: float, + point: BatteryPoint | None, + ) -> float: + if point is not None: + return point.current_a + return battery_power_w / max(1e-6, battery_voltage_v) + + @staticmethod + def _battery_point( + battery: BatterySpec, + current_a: float, + state: BatteryState | None, + ) -> BatteryPoint | None: + if hasattr(battery, "state_at_current"): + if state is None: + raise ValueError("battery_state is required for dynamic battery models") + return battery.state_at_current(state=state, current_a=current_a) + return None + + @staticmethod + def _battery_power( + battery: BatterySpec, + current_a: float, + motor_power_w: float, + system: SystemSpec, + ) -> float: + power_w = motor_power_w + (current_a**2) * system.resistance_ohm + discharge_efficiency = getattr(battery, "discharge_efficiency", 1.0) + return power_w / max(1e-6, discharge_efficiency) + + @staticmethod + def _battery_efficiency( + battery: BatterySpec, + point: BatteryPoint | None, + ) -> float: + if point is not None: + return point.efficiency + return float(getattr(battery, "discharge_efficiency", 1.0)) + @staticmethod def _estimate_j_max(prop_entry: PropellerEntry) -> float: j_maxes = [] diff --git a/scripts/generate_readme_plots.py b/scripts/generate_readme_plots.py index 1c1f237b..4208006b 100644 --- a/scripts/generate_readme_plots.py +++ b/scripts/generate_readme_plots.py @@ -1,38 +1,101 @@ +import sys +from pathlib import Path + import matplotlib -matplotlib.use('Agg') + +matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np -from pathlib import Path -import sys -# Add project root to sys.path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from pythrust.battery import BatteryState, FixedVoltageBattery, RateMapBattery from pythrust.propellers import PropellerDatabase -from pythrust.propulsion import ( - BatterySpec, - MotorSpec, - PropellerSpec, - SystemSpec, - PropulsionSolver, -) +from pythrust.propulsion import MotorSpec, PropellerSpec, PropulsionSolver, SystemSpec from pythrust.propulsion.autotune import ManufacturerTestPoint, PropulsionCalibrator -def generate_calibration_plot(db, prop_entry): + +IMAGE_DIR = Path("docs/images") +FIGSIZE = (9.6, 5.4) +DPI = 170 + +COLORS = { + "blue": "#2563eb", + "green": "#059669", + "orange": "#d97706", + "red": "#dc2626", + "purple": "#7c3aed", + "teal": "#0891b2", + "ink": "#111827", + "muted": "#6b7280", + "grid": "#d1d5db", + "panel": "#f8fafc", +} + + +def configure_style(): + plt.rcParams.update( + { + "figure.facecolor": "white", + "axes.facecolor": "white", + "axes.edgecolor": "#374151", + "axes.labelcolor": COLORS["ink"], + "axes.titlecolor": COLORS["ink"], + "xtick.color": COLORS["ink"], + "ytick.color": COLORS["ink"], + "grid.color": COLORS["grid"], + "grid.linewidth": 0.8, + "font.size": 10, + "axes.titlesize": 12, + "axes.labelsize": 10, + "legend.fontsize": 9, + "figure.titlesize": 15, + } + ) + + +def save_figure(fig, path): + fig.savefig(path, dpi=DPI) + plt.close(fig) + print(f"Saved: {path}") + + +def apply_two_panel_layout(fig): + fig.subplots_adjust( + left=0.09, + right=0.90, + bottom=0.17, + top=0.78, + wspace=0.48, + ) + + +def build_baseline_specs(): + motor = MotorSpec( + kv_rpm_per_v=860.0, + resistance_ohm=0.0258, + no_load_current_a=1.3, + current_max_a=65.0, + ) + system = SystemSpec(resistance_ohm=0.095) + propeller = PropellerSpec(diameter_m=0.3302) + return motor, system, propeller + + +def generate_calibration_plot(prop_entry): print("Generating Calibration Plot...") - motor = MotorSpec( kv_rpm_per_v=860.0, resistance_ohm=0.0258, no_load_current_a=1.3, current_max_a=65.0, ) - battery = BatterySpec(voltage_v=14.8) + battery = FixedVoltageBattery(voltage_v=14.8) system = SystemSpec(resistance_ohm=0.05) propeller = PropellerSpec(diameter_m=0.3302) - - RAW_TABLE = [ + + raw_table = [ {"rpm": 3897, "thrust_g": 500, "current_a": 3.9}, {"rpm": 4804, "thrust_g": 750, "current_a": 6.7}, {"rpm": 5421, "thrust_g": 1000, "current_a": 10.2}, @@ -45,14 +108,12 @@ def generate_calibration_plot(db, prop_entry): {"rpm": 8695, "thrust_g": 2750, "current_a": 47.5}, {"rpm": 9230, "thrust_g": 3350, "current_a": 63.2}, ] - points = [ ManufacturerTestPoint(rpm=r["rpm"], thrust_g=r["thrust_g"], current_a=r["current_a"]) - for r in RAW_TABLE + for r in raw_table ] - - calibrator = PropulsionCalibrator() - result = calibrator.calibrate( + + result = PropulsionCalibrator().calibrate( test_points=points, motor=motor, battery=battery, @@ -61,16 +122,13 @@ def generate_calibration_plot(db, prop_entry): prop_entry=prop_entry, ) fitted_system = result.to_system_spec() - solver = PropulsionSolver() - - # Generate smooth fitted predictions + model_rpms = np.linspace(3500, 9500, 100) model_thrusts_g = [] model_currents_a = [] - for rpm in model_rpms: - pt = solver._build_point( + point = solver._build_point( motor=motor, battery=battery, system=fitted_system, @@ -78,185 +136,297 @@ def generate_calibration_plot(db, prop_entry): prop_entry=prop_entry, rho=1.225, airspeed_mps=0.0, - rpm=rpm + rpm=rpm, ) - model_thrusts_g.append(pt.thrust_n * 1000.0 / 9.80665) - model_currents_a.append(pt.battery_power_w / battery.voltage_v) - - # Plotting - fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5)) - - # Left subplot: Thrust vs RPM - ax1.plot(model_rpms, model_thrusts_g, color='C0', label='Fitted Model') - ax1.scatter([p.rpm for p in points], [p.thrust_g for p in points], color='red', marker='o', label='Datasheet') - ax1.set_xlabel('RPM') - ax1.set_ylabel('Thrust (g)') - ax1.set_title('Thrust vs RPM') + model_thrusts_g.append(point.thrust_n * 1000.0 / 9.80665) + model_currents_a.append(point.battery_current_a) + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=FIGSIZE) + fig.suptitle("System Resistance Calibration", weight="bold") + fig.text( + 0.5, + 0.91, + f"Fitted wiring/ESC loss term: {result.system_resistance_ohm:.3f} ohm", + ha="center", + color=COLORS["muted"], + fontsize=10, + ) + + ax1.plot(model_rpms, model_thrusts_g, color=COLORS["blue"], linewidth=2.2, label="PyThrust fit") + ax1.scatter( + [p.rpm for p in points], + [p.thrust_g for p in points], + color=COLORS["red"], + s=32, + label="test data", + zorder=3, + ) + ax1.set_title("Thrust match") + ax1.set_xlabel("Shaft speed [RPM]") + ax1.set_ylabel("Thrust [g]") ax1.grid(True) - ax1.legend() - - # Right subplot: Current vs RPM - ax2.plot(model_rpms, model_currents_a, color='C1', label='Fitted Model') - ax2.scatter([p.rpm for p in points], [p.current_a for p in points], color='red', marker='o', label='Datasheet') - ax2.set_xlabel('RPM') - ax2.set_ylabel('Current (A)') - ax2.set_title('Current vs RPM') + ax1.legend(frameon=False) + + ax2.plot(model_rpms, model_currents_a, color=COLORS["orange"], linewidth=2.2, label="PyThrust fit") + ax2.scatter( + [p.rpm for p in points], + [p.current_a for p in points], + color=COLORS["red"], + s=32, + label="test data", + zorder=3, + ) + ax2.set_title("Battery current match") + ax2.set_xlabel("Shaft speed [RPM]") + ax2.set_ylabel("Battery current [A]") ax2.grid(True) - ax2.legend() - - fig.suptitle('Propulsion Calibration & Auto-Tuning Results', fontsize=14) - plt.tight_layout() - - out_path = Path("docs/images/calibration_results.png") - plt.savefig(out_path, dpi=150, bbox_inches='tight') - plt.close() - print(f"Saved: {out_path}") + ax2.legend(frameon=False) + + apply_two_panel_layout(fig) + save_figure(fig, IMAGE_DIR / "calibration_results.png") + def generate_propeller_plot(prop_entry): - print("Generating Propeller Coefficients Plot...") - - fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5)) - - # Plot Ct and Cp vs J for a few RPM bands + print("Generating Propeller Database Plot...") + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=FIGSIZE) + fig.suptitle("Empirical Propeller Database", weight="bold") + fig.text( + 0.5, + 0.91, + "APC 13x6.5E coefficient maps used by the operating-point solver", + ha="center", + color=COLORS["muted"], + fontsize=10, + ) + rpms = prop_entry.rpm_levels - # Pick a few evenly spaced RPM bands - selected_rpms = [rpms[i] for i in [0, len(rpms)//3, 2*len(rpms)//3, len(rpms)-1] if i < len(rpms)] - selected_rpms = sorted(list(set(selected_rpms))) - - for rpm in selected_rpms: + selected = [rpms[i] for i in [0, len(rpms) // 3, 2 * len(rpms) // 3, len(rpms) - 1] if i < len(rpms)] + selected = sorted(set(selected)) + palette = [COLORS["blue"], COLORS["green"], COLORS["orange"], COLORS["purple"]] + + for color, rpm in zip(palette, selected): points = prop_entry.data_by_rpm[rpm] j_vals = [p.j for p in points] - ct_vals = [p.ct for p in points] - cp_vals = [p.cp for p in points] - - ax1.plot(j_vals, ct_vals, label=f'{int(rpm)} RPM') - ax2.plot(j_vals, cp_vals, label=f'{int(rpm)} RPM') - - ax1.set_xlabel('Advance Ratio (J)') - ax1.set_ylabel('Thrust Coefficient (Ct)') - ax1.set_title('Ct vs J') + ax1.plot(j_vals, [p.ct for p in points], color=color, linewidth=2.0, label=f"{int(rpm)} RPM") + ax2.plot(j_vals, [p.cp for p in points], color=color, linewidth=2.0, label=f"{int(rpm)} RPM") + + ax1.set_title("Thrust coefficient") + ax1.set_xlabel("Advance ratio J") + ax1.set_ylabel("Ct") + ax1.grid(True) + ax1.legend(frameon=False) + + ax2.set_title("Power coefficient") + ax2.set_xlabel("Advance ratio J") + ax2.set_ylabel("Cp") + ax2.grid(True) + ax2.legend(frameon=False) + + apply_two_panel_layout(fig) + save_figure(fig, IMAGE_DIR / "propeller_coefficients.png") + + +def generate_battery_mission_plot(prop_entry): + print("Generating Rate-Map Battery Mission Plot...") + battery = RateMapBattery.from_json("data/batteries/example_liion_cell.json", series=4, parallel=2) + motor, system, propeller = build_baseline_specs() + solver = PropulsionSolver() + state = BatteryState(soc=0.95) + + segments = [ + {"name": "takeoff", "duration_s": 45.0, "throttle": 0.70, "airspeed_mps": 0.0}, + {"name": "climb", "duration_s": 90.0, "throttle": 0.65, "airspeed_mps": 5.0}, + {"name": "cruise", "duration_s": 180.0, "throttle": 0.50, "airspeed_mps": 10.0}, + {"name": "return", "duration_s": 120.0, "throttle": 0.45, "airspeed_mps": 8.0}, + ] + + times_min = [0.0] + soc = [state.soc] + voltages = [] + currents = [] + c_rates = [] + labels = [] + elapsed_s = 0.0 + + for segment in segments: + point = solver.solve_operating_point( + motor=motor, + battery=battery, + battery_state=state, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=segment["airspeed_mps"], + throttle=segment["throttle"], + ) + elapsed_s += segment["duration_s"] + state = battery.step_current(state=state, current_a=point.battery_current_a, dt_s=segment["duration_s"]) + times_min.append(elapsed_s / 60.0) + soc.append(state.soc) + voltages.append(point.battery_voltage_v) + currents.append(point.battery_current_a) + c_rates.append(point.battery_c_rate) + labels.append(segment["name"]) + + centers = [(times_min[i] + times_min[i + 1]) / 2 for i in range(len(segments))] + widths = [times_min[i + 1] - times_min[i] for i in range(len(segments))] + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=FIGSIZE) + fig.suptitle("Rate-Map Battery Mission Example", weight="bold") + fig.text( + 0.5, + 0.91, + "A short segment schedule updates battery state from solved pack current", + ha="center", + color=COLORS["muted"], + fontsize=10, + ) + + ax1.plot(times_min, soc, color=COLORS["green"], marker="o", linewidth=2.3, label="SoC") + ax1.set_title("State of charge") + ax1.set_xlabel("Elapsed time [min]") + ax1.set_ylabel("SoC") + ax1.set_ylim(0.75, 0.98) ax1.grid(True) - ax1.legend() - - ax2.set_xlabel('Advance Ratio (J)') - ax2.set_ylabel('Power Coefficient (Cp)') - ax2.set_title('Cp vs J') + + ax2.plot( + centers, + currents, + color=COLORS["orange"], + marker="o", + linewidth=2.4, + label="current", + ) + ax2.set_title("Solved pack current") + ax2.set_xlabel("Mission segment") + ax2.set_ylabel("Battery current [A]") + ax2.set_xticks(centers) + ax2.set_xticklabels(labels, rotation=20, ha="right") ax2.grid(True) - ax2.legend() - - fig.suptitle('Propeller Aerodynamic Coefficients (APC 13x6.5E)', fontsize=14) - plt.tight_layout() - - out_path = Path("docs/images/propeller_coefficients.png") - plt.savefig(out_path, dpi=150, bbox_inches='tight') - plt.close() - print(f"Saved: {out_path}") + ax2.set_ylim(0.0, max(currents) * 1.18) + for x, current, segment in zip(centers, currents, segments): + ax2.text( + x, + current + max(currents) * 0.035, + f"{segment['throttle'] * 100:.0f}%", + ha="center", + va="bottom", + color=COLORS["muted"], + fontsize=9, + ) + + apply_two_panel_layout(fig) + save_figure(fig, IMAGE_DIR / "rate_map_battery_mission.png") + def generate_heatmap_plot(db): print("Generating Efficiency Heatmap...") - - # Set up a grid of Motor Kv vs Propeller Diameter kv_grid = np.linspace(400, 1200, 25) - dia_grid = np.linspace(10, 18, 25) # in inches - - Z_eff = np.zeros((len(dia_grid), len(kv_grid))) - - battery = BatterySpec(voltage_v=14.8) + dia_grid = np.linspace(10, 18, 25) + efficiency = np.zeros((len(dia_grid), len(kv_grid))) + + battery = FixedVoltageBattery(voltage_v=14.8) system = SystemSpec(resistance_ohm=0.05) solver = PropulsionSolver() - - # Target 500 grams of thrust (4.903 N) target_thrust_n = 4.903 - + def find_hover_throttle(motor, propeller, prop_entry): - def residual(throttle): - pt = solver.solve_operating_point( - motor=motor, battery=battery, system=system, - propeller=propeller, prop_entry=prop_entry, - rho=1.225, airspeed_mps=0.0, throttle=throttle - ) - if not pt.is_feasible: - # Retain directionality - return pt.thrust_n - target_thrust_n - return pt.thrust_n - target_thrust_n - - # Binary search for throttle low, high = 0.1, 0.99 + mid = 0.5 for _ in range(15): mid = (low + high) / 2.0 - res = residual(mid) - if abs(res) < 1e-3: + point = solver.solve_operating_point( + motor=motor, + battery=battery, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=0.0, + throttle=mid, + ) + residual = point.thrust_n - target_thrust_n + if abs(residual) < 1e-3: break - if res < 0: + if residual < 0: low = mid else: high = mid - - pt = solver.solve_operating_point( - motor=motor, battery=battery, system=system, - propeller=propeller, prop_entry=prop_entry, - rho=1.225, airspeed_mps=0.0, throttle=mid + + point = solver.solve_operating_point( + motor=motor, + battery=battery, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=0.0, + throttle=mid, ) - return pt if (pt.is_feasible and abs(pt.thrust_n - target_thrust_n) < 0.1) else None + return point if point.is_feasible and abs(point.thrust_n - target_thrust_n) < 0.1 else None for i, dia_in in enumerate(dia_grid): - # Retrieve the propeller entry closest to this size prop_entry = db.find_by_size(dia_in, pitch_in=dia_in * 0.5, blade_count=2, tolerance=2.0) if prop_entry is None: - # Fallback to standard 13x6.5E prop_entry = db.get("APC_13x6.5E") - propeller = PropellerSpec(diameter_m=dia_in * 0.0254) - + for j, kv in enumerate(kv_grid): motor = MotorSpec( kv_rpm_per_v=kv, resistance_ohm=0.0258, no_load_current_a=1.3, - current_max_a=65.0 + current_max_a=65.0, ) - - pt = find_hover_throttle(motor, propeller, prop_entry) - if pt is not None and pt.battery_power_w > 0: - thrust_g = pt.thrust_n * 1000.0 / 9.80665 - Z_eff[i, j] = thrust_g / pt.battery_power_w # g/W + point = find_hover_throttle(motor, propeller, prop_entry) + if point is not None and point.battery_power_w > 0: + thrust_g = point.thrust_n * 1000.0 / 9.80665 + efficiency[i, j] = thrust_g / point.battery_power_w else: - Z_eff[i, j] = np.nan - - # Plotting - plt.figure(figsize=(7.5, 5.5)) - cp = plt.contourf(kv_grid, dia_grid, Z_eff, levels=15, cmap='viridis') - cbar = plt.colorbar(cp) - cbar.set_label('Hover Efficiency (g/W)', rotation=270, labelpad=15) - - plt.xlabel('Motor Kv (RPM/V)') - plt.ylabel('Propeller Diameter (in)') - plt.title('Propulsion Hover Efficiency Map (at 500g Thrust)') - plt.grid(True, linestyle=':', alpha=0.6) - - out_path = Path("docs/images/efficiency_heatmap.png") - plt.savefig(out_path, dpi=150, bbox_inches='tight') - plt.close() - print(f"Saved: {out_path}") + efficiency[i, j] = np.nan + + fig = plt.figure(figsize=FIGSIZE) + ax = fig.add_axes([0.08, 0.13, 0.76, 0.70]) + contour = ax.contourf(kv_grid, dia_grid, efficiency, levels=16, cmap="viridis") + cax = fig.add_axes([0.87, 0.13, 0.025, 0.70]) + cbar = fig.colorbar(contour, cax=cax) + cbar.set_label("Hover efficiency [g/W]", rotation=270, labelpad=15) + ax.set_xlabel("Motor Kv [RPM/V]") + ax.set_ylabel("Propeller diameter [in]") + fig.suptitle("Hover Efficiency Design Map", fontweight="bold") + fig.text( + 0.5, + 0.91, + "Target thrust: 500 g", + ha="center", + color=COLORS["muted"], + fontsize=10, + ) + ax.grid(True, linestyle=":", alpha=0.6) + save_figure(fig, IMAGE_DIR / "efficiency_heatmap.png") + def main(): - # Load databases + configure_style() + db_prop = PropellerDatabase() if not db_prop.load(Path("data/propellers/apc_202602"), strict=False): print("Error: Could not load propeller database.") sys.exit(1) - + prop_entry = db_prop.get("APC_13x6.5E") if prop_entry is None: print("Error: Propeller APC_13x6.5E not found.") sys.exit(1) - - # Ensure docs/images directory exists - Path("docs/images").mkdir(parents=True, exist_ok=True) - - generate_calibration_plot(db_prop, prop_entry) + + IMAGE_DIR.mkdir(parents=True, exist_ok=True) + + generate_battery_mission_plot(prop_entry) + generate_calibration_plot(prop_entry) generate_propeller_plot(prop_entry) generate_heatmap_plot(db_prop) -if __name__ == '__main__': + +if __name__ == "__main__": main() diff --git a/tests/test_battery.py b/tests/test_battery.py new file mode 100644 index 00000000..51f382c6 --- /dev/null +++ b/tests/test_battery.py @@ -0,0 +1,153 @@ +import json +import math +from pathlib import Path + +import pytest + +from pythrust.battery import BatteryState, FixedVoltageBattery, RateMapBattery +from pythrust.propulsion.models import BatterySpec + + +PROJECT_ROOT = Path(__file__).resolve().parents[1] + + +@pytest.fixture +def rate_map_battery(): + return RateMapBattery( + name="Example pack", + capacity_ah=4.2, + cutoff_voltage_v=2.5, + charge_voltage_v=4.2, + max_current_a=20.0, + series=4, + parallel=2, + dod=[0.0, 0.5, 1.0], + ocv_v=[4.2, 3.8, 3.2], + resistance_ohm=[0.02, 0.03, 0.05], + ) + + +def test_fixed_voltage_battery_replaces_battery_spec(): + battery = FixedVoltageBattery(voltage_v=14.8, discharge_efficiency=0.95) + + assert BatterySpec is FixedVoltageBattery + assert battery.terminal_voltage(current_a=10.0) == 14.8 + assert math.isclose(battery.terminal_power(current_a=10.0), 14.8 * 10.0 / 0.95) + + +def test_battery_state_dod_round_trip(): + state = BatteryState.from_dod(0.25) + + assert state.soc == 0.75 + assert state.dod == 0.25 + + +@pytest.mark.parametrize("soc", [math.nan, math.inf, -math.inf]) +def test_battery_state_rejects_non_finite_soc(soc): + with pytest.raises(ValueError, match="soc"): + BatteryState(soc=soc) + + +@pytest.mark.parametrize("dod", [math.nan, math.inf, -math.inf]) +def test_battery_state_rejects_non_finite_dod(dod): + with pytest.raises(ValueError, match="dod"): + BatteryState.from_dod(dod) + + +def test_rate_map_state_at_current(rate_map_battery): + state = BatteryState(soc=0.5) + point = rate_map_battery.state_at_current(state=state, current_a=8.0) + + assert point.is_feasible is True + assert math.isclose(point.cell_current_a, 4.0) + assert math.isclose(point.cell_voltage_v, 3.8 - 0.03 * 4.0) + assert math.isclose(point.terminal_voltage_v, (3.8 - 0.03 * 4.0) * 4) + assert math.isclose(point.power_w, point.terminal_voltage_v * 8.0) + assert math.isclose(point.c_rate, 4.0 / 4.2) + + +def test_rate_map_state_at_power(rate_map_battery): + state = BatteryState(soc=0.5) + point = rate_map_battery.state_at_power(state=state, power_w=100.0) + + ocv = 3.8 + resistance = 0.03 + cell_power = 100.0 / (4 * 2) + expected_cell_current = (ocv - math.sqrt(ocv**2 - 4.0 * resistance * cell_power)) / ( + 2.0 * resistance + ) + + assert point.is_feasible is True + assert math.isclose(point.cell_current_a, expected_cell_current) + assert math.isclose(point.power_w, 100.0, rel_tol=1e-12) + + +def test_rate_map_reports_infeasible_power(rate_map_battery): + state = BatteryState(soc=0.5) + point = rate_map_battery.state_at_power(state=state, power_w=2000.0) + + assert point.is_feasible is False + assert point.infeasible_reason == "power_limit" + + +def test_rate_map_preserves_first_infeasible_reason(rate_map_battery): + state = BatteryState(soc=0.5) + point = rate_map_battery.state_at_current(state=state, current_a=1000.0) + + assert point.is_feasible is False + assert point.cell_current_a > rate_map_battery.max_current_a + assert point.cell_voltage_v < rate_map_battery.cutoff_voltage_v + assert point.infeasible_reason == "current_limit" + + +def test_rate_map_step_current(rate_map_battery): + state = BatteryState(soc=1.0) + next_state = rate_map_battery.step_current(state=state, current_a=8.4, dt_s=1800.0) + + assert math.isclose(next_state.dod, 0.5) + assert math.isclose(next_state.soc, 0.5) + + +def test_rate_map_loads_from_json(tmp_path): + path = tmp_path / "battery.json" + path.write_text( + json.dumps( + { + "name": "JSON pack", + "source": "test", + "cell": { + "capacity_ah": 4.2, + "cutoff_voltage_v": 2.5, + "charge_voltage_v": 4.2, + "max_current_a": 20.0, + }, + "curves": { + "dod": [0.0, 0.5, 1.0], + "ocv_v": [4.2, 3.8, 3.2], + "resistance_ohm": [0.02, 0.03, 0.05], + }, + } + ), + encoding="utf-8", + ) + + battery = RateMapBattery.from_json(path, series=4, parallel=2) + + assert battery.name == "JSON pack" + assert battery.source == "test" + assert battery.series == 4 + assert battery.parallel == 2 + + +def test_rate_map_loads_example_cell_dataset(): + path = PROJECT_ROOT / "data" / "batteries" / "example_liion_cell.json" + + battery = RateMapBattery.from_json(path, series=4, parallel=2) + point = battery.state_at_current(state=BatteryState(soc=0.5), current_a=8.4) + + assert battery.name == "Example Li-ion Cell" + assert battery.series == 4 + assert battery.parallel == 2 + assert battery.capacity_ah == 4.2 + assert point.is_feasible is True + assert math.isclose(point.cell_current_a, 4.2) diff --git a/tests/test_models.py b/tests/test_models.py index 5680400d..f07b3d91 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -1,12 +1,10 @@ import math import pytest -from pythrust.propulsion.models import ( - MotorSpec, - BatterySpec, - SystemSpec, - PropellerSpec, - OperatingPoint -) +from pythrust.propulsion.models.battery import BatterySpec +from pythrust.propulsion.models.motor import MotorSpec +from pythrust.propulsion.models.operating_point import OperatingPoint +from pythrust.propulsion.models.propeller import PropellerSpec +from pythrust.propulsion.models.system import SystemSpec def test_motor_spec_no_load_current(): @@ -105,7 +103,11 @@ def test_operating_point(): is_feasible=True, propeller_efficiency=0.65, motor_efficiency=0.83, - system_efficiency=0.5395 + system_efficiency=0.5395, + battery_voltage_v=11.1, + battery_current_a=27.9, + battery_c_rate=6.64, + battery_efficiency=0.98, ) assert op.rpm == 8000.0 assert op.is_feasible is True @@ -113,6 +115,10 @@ def test_operating_point(): assert op.propeller_efficiency == 0.65 assert op.motor_efficiency == 0.83 assert op.system_efficiency == 0.5395 + assert op.battery_voltage_v == 11.1 + assert op.battery_current_a == 27.9 + assert op.battery_c_rate == 6.64 + assert op.battery_efficiency == 0.98 # Test default values op_default = OperatingPoint( @@ -132,4 +138,7 @@ def test_operating_point(): assert op_default.propeller_efficiency == 0.0 assert op_default.motor_efficiency == 0.0 assert op_default.system_efficiency == 0.0 - + assert op_default.battery_voltage_v == 0.0 + assert op_default.battery_current_a == 0.0 + assert op_default.battery_c_rate == 0.0 + assert op_default.battery_efficiency == 1.0 diff --git a/tests/test_solver.py b/tests/test_solver.py index 429c73e5..4670a214 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -1,5 +1,6 @@ import math import pytest +from pythrust.battery import BatteryState, RateMapBattery from pythrust.propellers.database import ( PropellerMetadata, PropellerDataPoint, @@ -75,6 +76,10 @@ def test_solver_zero_throttle(test_setup): assert op.rpm == 0.0 assert op.is_feasible is False assert op.infeasible_reason == "throttle<=0" + assert op.battery_voltage_v == battery.voltage_v + assert op.battery_current_a == 0.0 + assert op.battery_c_rate == 0.0 + assert op.battery_efficiency == battery.discharge_efficiency def test_solver_successful_solve(test_setup): @@ -109,6 +114,137 @@ def test_solver_successful_solve(test_setup): v_motor = op.motor_voltage_v i_sys = op.motor_current_a assert math.isclose(v_motor + i_sys * system.resistance_ohm, v_applied, rel_tol=1e-3) + assert math.isclose(op.battery_voltage_v, battery.voltage_v) + assert math.isclose(op.battery_current_a, op.battery_power_w / battery.voltage_v) + assert op.battery_c_rate == 0.0 + assert op.battery_efficiency == battery.discharge_efficiency + + +def test_solver_rate_map_matches_fixed_voltage_when_voltage_is_flat(test_setup): + motor, battery, system, propeller, prop_entry = test_setup + solver = PropulsionSolver() + state = BatteryState(soc=0.5) + rate_map = RateMapBattery( + name="flat voltage", + capacity_ah=4.2, + cutoff_voltage_v=2.5, + charge_voltage_v=4.2, + max_current_a=1000.0, + series=3, + parallel=1, + dod=[0.0, 1.0], + ocv_v=[battery.voltage_v / 3.0, battery.voltage_v / 3.0], + resistance_ohm=[1e-12, 1e-12], + ) + + fixed_op = 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, + ) + rate_map_op = solver.solve_operating_point( + motor=motor, + battery=rate_map, + battery_state=state, + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=0.0, + throttle=0.8, + ) + + assert rate_map_op.is_feasible is True + assert math.isclose(rate_map_op.rpm, fixed_op.rpm, rel_tol=1e-6) + assert math.isclose(rate_map_op.battery_voltage_v, battery.voltage_v, rel_tol=1e-6) + assert math.isclose(rate_map_op.battery_current_a, 0.8 * rate_map_op.motor_current_a, rel_tol=1e-6) + assert rate_map_op.battery_c_rate > 0.0 + assert math.isclose(rate_map_op.battery_efficiency, 1.0, rel_tol=1e-9) + + +def test_solver_rate_map_voltage_changes_with_state_of_charge(test_setup): + motor, _, system, propeller, prop_entry = test_setup + solver = PropulsionSolver() + battery = RateMapBattery( + name="soc-sensitive voltage", + capacity_ah=4.2, + cutoff_voltage_v=2.5, + charge_voltage_v=4.2, + max_current_a=1000.0, + series=3, + parallel=2, + dod=[0.0, 1.0], + ocv_v=[4.1, 3.4], + resistance_ohm=[0.01, 0.03], + ) + + high_soc = solver.solve_operating_point( + motor=motor, + battery=battery, + battery_state=BatteryState(soc=1.0), + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=0.0, + throttle=0.8, + ) + low_soc = solver.solve_operating_point( + motor=motor, + battery=battery, + battery_state=BatteryState(soc=0.2), + system=system, + propeller=propeller, + prop_entry=prop_entry, + rho=1.225, + airspeed_mps=0.0, + throttle=0.8, + ) + + assert high_soc.is_feasible is True + assert low_soc.is_feasible is True + assert high_soc.rpm > low_soc.rpm + assert high_soc.battery_voltage_v > low_soc.battery_voltage_v + assert high_soc.battery_current_a > 0.0 + assert low_soc.battery_current_a > 0.0 + assert high_soc.battery_c_rate > 0.0 + assert low_soc.battery_c_rate > 0.0 + assert 0.0 < high_soc.battery_efficiency <= 1.0 + assert 0.0 < low_soc.battery_efficiency <= 1.0 + + +def test_solver_rate_map_requires_battery_state(test_setup): + motor, _, system, propeller, prop_entry = test_setup + solver = PropulsionSolver() + battery = RateMapBattery( + name="requires state", + capacity_ah=4.2, + cutoff_voltage_v=2.5, + charge_voltage_v=4.2, + max_current_a=1000.0, + series=3, + parallel=1, + dod=[0.0, 1.0], + ocv_v=[4.1, 3.4], + resistance_ohm=[0.01, 0.03], + ) + + with pytest.raises(ValueError, match="battery_state"): + 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, + ) def test_solver_current_limit_exceeded(test_setup): @@ -263,3 +399,45 @@ def test_solver_invalid_efficiency(test_setup): assert op.infeasible_reason == "invalid_efficiency" +def test_solver_preserves_first_infeasible_reason(test_setup): + motor, battery, system, propeller, _ = test_setup + low_limit_motor = MotorSpec( + kv_rpm_per_v=motor.kv_rpm_per_v, + resistance_ohm=motor.resistance_ohm, + no_load_current_a=motor.no_load_current_a, + current_max_a=0.01, + ) + 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 = { + 6000.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() + rpm = 5.0 / (0.2 * propeller.diameter_m) * 60.0 + + op = solver._build_point( + motor=low_limit_motor, + battery=battery, + system=system, + propeller=propeller, + prop_entry=bad_entry, + rho=1.225, + airspeed_mps=5.0, + rpm=rpm, + throttle=0.8, + ) + + assert op.is_feasible is False + assert op.infeasible_reason == "current_limit"