Improved mass iteration schemes

docs/src/trajectories/trajectories.md

@@ -19,6 +19,7 @@ vertical, speed) trajectories.
 :maxdepth: 1
 
 trajectory_data.md
+trajectory_mass_iteration.md
 trajectory_builders.md
 trajectory_stores.md
 ```

docs/src/trajectories/trajectory_mass_iteration.md

@@ -0,0 +1,65 @@
+# Iterating Fuel Mass During Trajectory Simulation
+
+The v0 implementation of mass iteration in the new AEIC is based on a simple naive correction with a set maximum number of iteration. In this case, the mass residual is specified as the actual fuel burned over the predicted total fuel mass:
+
+$$
+\delta = \frac{\text{Actual Fuel Burn} - M_f}{M_f}
+$$
+
+This correction can be improved using an approximation from the Breguet range equation. Assuming the aerodynamic and engine properties are negligibly changed by the mass variation, the original `_fly_iteration` call used $\delta$ percent more fuel than was loaded and flew a range, $R$:
+
+$$
+R \propto -c^*\ln\bigg(1 - \frac{M_f}{M_0}(1+\delta)\bigg)
+$$
+
+We call the initial iteration's fuel-to-initial mass fraction $\lambda$ such that:
+
+$$
+R \propto -c^*\ln\bigg(1 - \lambda(1+\delta)\bigg)
+$$
+
+We want to fly the same mission (same range) with some additional amount of fuel relative to the initial prediction, $\Delta = \text{Additional Fuel}/M_f$. This would correspond to:
+
+$$
+R \propto -c^*\ln\bigg(1 - \frac{1 + \Delta}{1/\lambda + \Delta}\bigg)
+$$
+
+Again, we assume the aero and engine properties remain constant such that $c^*\equiv const$. Using this, we set the ranges equal and get:
+
+$$
+1-\lambda(1+\delta) = 1 - \frac{1 +\Delta}{1/\lambda + \Delta}
+$$
+
+$$
+\frac{1 + \Delta}{1/\lambda + \Delta} = \lambda(1+\delta)
+$$
+
+$$
+1 + \Delta = \bigg(\frac{1}{\lambda} + \Delta\bigg)\lambda(1 + \delta)
+$$
+
+$$
+1 + \Delta = (1 + \lambda\Delta)(1 + \delta)
+$$
+
+$$
+1 + \Delta = 1 + \lambda\Delta + \delta + \delta\lambda\Delta
+$$
+
+$$
+\Delta(1 - \lambda - \delta\lambda) = \delta
+$$
+
+$$
+\Delta = \frac{\delta}{1 - \lambda(1 + \delta)}
+$$
+
+In the naive approach, we simply set $\Delta = \delta$. The difference between these methods is shown in the figure below:
+
+```{image} ../../_static/Naive-vs-rangecorr.png
+:align: center
+```
+
+We see that for cases where additional fuel is needed to complete the mission ($\delta > 0$), we must add more fuel than the naive approach suggests. Conversely, for missions with too much initial fuel ($\delta<0$), we can remove more fuel than predicted by the naive method.
+
+However, when testing these implementations in `notebooks/mass_iteration.ipynb`, the range-corrected method was found to be overshooting the correction, leading to increased iterations. To address this, a simple arithmetic average of the naive and range-corrected methods is preferred.

src/AEIC/trajectories/builders/__init__.py

@@ -1,6 +1,6 @@
 from .adjustable_legacy import AdjustableLegacyBuilder
 from .ads_b import ADSBBuilder, ADSBOptions
-from .base import Builder, Context, Options
+from .base import Builder, Context, MassIterMethod, Options
 from .dymos import DymosBuilder, DymosOptions
 from .legacy import LegacyBuilder, LegacyOptions
 from .tasopt import TASOPTBuilder, TASOPTOptions
@@ -9,6 +9,7 @@
     'Builder',
     'Context',
     'Options',
+    'MassIterMethod',
     'TASOPTBuilder',
     'TASOPTOptions',
     'ADSBBuilder',

src/AEIC/trajectories/builders/base.py

@@ -1,6 +1,7 @@
 import logging
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
+from enum import Enum, auto
 
 from AEIC.missions import Mission
 from AEIC.performance.model_selector import PerformanceModelSelector
@@ -14,6 +15,25 @@
 logger = logging.getLogger(__name__)
 
 
+class MassIterMethod(Enum):
+    """Type of mass iteration update scheme to use.
+
+    - NAIVE - If the mission ends with $dM$ more fuel burned than we initially
+        loaded, adjust the starting mass and fuel mass by $+dM$ and try flying the
+        mission again.
+    - RANGE_CORR - Using the range equation, account for the additional fuel
+        required to carry the $+dM$ extra while meeting the same range.
+    - AVERAGED - A simple arithmetic average of Naive and Range-Corrected
+        approaches. This is due to the Range-Corrected method tending to overestimate
+        the amount of fuel to be added/subtracted, leading to increased iterations vs.
+        the naive approach.
+    """
+
+    NAIVE = auto()
+    RANGE_CORR = auto()
+    AVERAGED = auto()
+
+
 @dataclass
 class Options:
     """Common options for trajectory builders."""
@@ -30,11 +50,14 @@ class Options:
     use_weather: bool = False
     """Whether to use wind data for ground-speed calculations."""
 
-    max_mass_iters: int = 5
-    """Maximum number of mass iterations (if used). Defaults to 5."""
+    iter_method: MassIterMethod = MassIterMethod.NAIVE
+    """Whether to use the naive or range-corrected mass iteration scheme."""
+
+    max_mass_iters: int = 10
+    """Maximum number of mass iterations (if used). Defaults to 10."""
 
-    mass_iter_reltol: float = 1e-2
-    """Desired relative tolerance for mass iteration. Defaults to 1e-2."""
+    mass_iter_reltol: float = 1e-6
+    """Desired relative tolerance for mass iteration. Defaults to 1e-6."""
 
 
 @dataclass
@@ -206,7 +229,10 @@ def fly(
 
             if self.options.iterate_mass:
                 # Iterate on starting mass to minimize mass residual.
-                traj = self._iterate_mass()
+                traj, iters = self._iterate_mass()
+
+                # Store the number of iterations required to converge.
+                self.iters = iters
             else:
                 # Otherwise, just fly a single iteration with the given starting
                 # mass.
@@ -250,9 +276,9 @@ def _iterate_mass(self) -> Trajectory:
             if abs(mass_res) < self.options.mass_iter_reltol:
                 mass_converged = True
             else:
-                # Perform a "dumb" correction of the starting mass.
-                self.starting_mass -= mass_res * self.total_fuel_mass
-                self.total_fuel_mass -= mass_res * self.total_fuel_mass
+                # Perform a correction of the starting mass.
+                self.starting_mass += mass_res * self.total_fuel_mass
+                self.total_fuel_mass += mass_res * self.total_fuel_mass
 
                 traj, mass_res = self._fly_iteration()
                 iter += 1
@@ -263,7 +289,7 @@ def _iterate_mass(self) -> Trajectory:
                 f"{mass_res:.2e} > {self.options.mass_iter_reltol:.2e}"
             )
 
-        return traj
+        return traj, iter
 
     def _start_point(self, traj: Trajectory) -> Container:
         # Set initial values, taking initial position and azimuth from ground
@@ -314,7 +340,18 @@ def _fly_iteration(self) -> tuple[Trajectory, float]:
 
         # Calculate weight residual normalized by total_fuel_mass.
         fuelBurned = self.starting_mass - traj.aircraft_mass[-1]
-        mass_residual = (self.total_fuel_mass - fuelBurned) / self.total_fuel_mass
+        delta = (fuelBurned - self.total_fuel_mass) / self.total_fuel_mass
+        lamda = self.total_fuel_mass / self.starting_mass
+
+        # Determine the update step to take for mass iteration based on the
+        # naive, range-correected, or averaged approach.
+        match self.options.iter_method:
+            case MassIterMethod.NAIVE:
+                mass_residual = delta
+            case MassIterMethod.RANGE_CORR:
+                mass_residual = delta / (1 - lamda * (1 + delta))
+            case MassIterMethod.AVERAGED:
+                mass_residual = 0.5 * ((delta / (1 - lamda * (1 + delta))) + delta)
 
         return traj, mass_residual