contrailopt.optimize

Trajectory optimization with Poll-Schumann Aircraft Performance and climate term.

Functions

cruise_flight_levels(origin_icao, dest_icao)

Determine the candidate cruise flight levels for a given origin-destination pair.

estimate_flight_hours(origin, dest[, ...])

Estimate upper-bound flight duration in hours.

solve_dag(dag, amass_init, fl_choices, ...)

Solve shortest-path DP on the topo-sorted DAG, tracking mass exactly.

solve_track(dag, profile, amass_init, ...)

Solve the vertical profile along a fixed track, choosing the flight level and Mach number.

Classes

DAGResult(*, state, amass_init, trip_fuel, ...)

The output of Optimizer.solve().

DAGState(*, best_cost, best_mass, best_time, ...)

Working arrays for the wavefront dynamic program.

Optimizer(origin_icao, dest_icao, ...[, ...])

Trajectory optimizer for a single origin-destination pair based on the PS model.

class contrailopt.optimize.DAGState(*, best_cost: NDArray[float32], best_mass: NDArray[float32], best_time: NDArray[float32], best_mach: NDArray[float32], best_climb_dist: NDArray[float32], best_climb_time: NDArray[float32], best_prev_h: NDArray[int64], best_prev_fi: NDArray[int64])[source]

Bases: object

Working arrays for the wavefront dynamic program.

classmethod initialize(n_h: int, n_cols: int) Self[source]

Initialize arrays to default values for a DAG search.

class contrailopt.optimize.DAGResult(*, state: DAGState, amass_init: float, trip_fuel: float, payload: float, reserve_fuel: float, landing_mass: float)[source]

Bases: object

The output of Optimizer.solve().

contrailopt.optimize.solve_dag(dag: HorizontalDAG, amass_init: float, fl_choices: NDArray[float32], mach_choices: NDArray[float32], atyp: PSAircraftEngineParams, cost_index: float, eef_cost_factor: float, origin_elev_ft: float, dest_elev_ft: float, takeoff_time: Timestamp, met_lookup: EdgeMetLookup | None, allow_cooling_credit: bool, step_penalty_kg: float, on_wavefront: Callable[[NDArray[int64], DAGState], None] | None = None) DAGState[source]

Solve shortest-path DP on the topo-sorted DAG, tracking mass exactly.

contrailopt.optimize.solve_track(dag: HorizontalDAG | Track, profile: Dataset, amass_init: float, fl_choices: NDArray[float32], mach_choices: NDArray[float32], atyp: PSAircraftEngineParams, cost_index: float, eef_cost_factor: float, origin_elev_ft: float, dest_elev_ft: float, allow_cooling_credit: bool, step_penalty_kg: float) DAGState[source]

Solve the vertical profile along a fixed track, choosing the flight level and Mach number.

Unlike solve_dag(), there is no precomputed edge set. From a given state, each action determines how far the aircraft advances: a level cruise (target FL equal to or below the current one) moves one waypoint; a climb moves as far as the climb model’s own distance requires, then cruises the remainder of that segment to the next waypoint.

contrailopt.optimize.estimate_flight_hours(origin: AirportCoords, dest: AirportCoords, mach_number: float = 0.75, max_headwind: float = 40.0) int[source]

Estimate upper-bound flight duration in hours.

Computes the worst-case flight time assuming the aircraft flies at mach_number at FL400 with a sustained headwind of max_headwind.

Parameters:
  • origin (AirportCoords) – Origin airport.

  • dest (AirportCoords) – Destination airport.

  • mach_number (float, default 0.75) – Cruise Mach number. Use the slowest aircraft’s design Mach minus a margin.

  • max_headwind (float, default 40.0) – Assumed maximum sustained headwind in m/s.

Returns:

Ceiling of estimated flight time in hours.

Return type:

int

contrailopt.optimize.cruise_flight_levels(origin_icao: str | AirportCoords, dest_icao: str | AirportCoords) NDArray[float32][source]

Determine the candidate cruise flight levels for a given origin-destination pair.

This function applies the common eastbound/westbound FL rules of even FLs for westbound flights and odd FLs for eastbound flights. There is not per-aircraft-type ceiling applied (this could be added if needed).

Parameters:
  • origin_icao (str | AirportCoords) – ICAO code for the origin airport (e.g. "KLAX") or pre-fetched coordinates.

  • dest_icao (str | AirportCoords) – ICAO code for the destination airport or pre-fetched coordinates.

Returns:

Array of candidate cruise flight levels in feet (e.g. [29000., 31000., ..., 41000.]).

Return type:

numpy.ndarray

class contrailopt.optimize.Optimizer(origin_icao: str | AirportCoords, dest_icao: str | AirportCoords, aircraft_type: str, takeoff_time: Timestamp, *, met: MetDataset | Dataset | None = None, eef: DataArray | MetDataArray | None = None, dag: HorizontalDAG | Track | None = None, cost_index: float = 60.0, dollar_tonne_co2e: float = 0.0, dollar_kg_fuel: float = 1.0, step_penalty_kg: float = 0.0, met_spacing_m: float = 25000.0, flight_hours: int | None = None, allow_cooling_credit: bool = False, avoidance_regions: list[list[tuple[float, float]]] | None = None, fl_choices: NDArray[float32] | None = None, **kwargs: Any)[source]

Bases: object

Trajectory optimizer for a single origin-destination pair based on the PS model.

The optimizer builds a horizontal directed acyclic graph between two airports, optionally interpolates met data onto edge sample points, and solves for the minimum-cost path across valid flight levels and Mach numbers.

Parameters:
  • origin_icao (str | AirportCoords) – ICAO code for the origin airport (e.g. "KLAX"), or an AirportCoords instance.

  • dest_icao (str | AirportCoords) – ICAO code for the destination airport, or an AirportCoords instance.

  • aircraft_type (str) – Aircraft type key in the PS model parameter table (e.g. "A320").

  • takeoff_time (pandas.Timestamp) – Departure time, used for met interpolation.

  • met (MetDataset | xr.Dataset | None, default None) – Gridded met data with air_temperature, eastward_wind, and northward_wind. If None, cruise performance uses ISA temperatures and zero wind.

  • eef (xr.DataArray | MetDataArray | None, default None) – Optional eef_per_m DataArray on its own lon/lat grid. If provided, EEF is interpolated onto sample points independently from the weather grid, avoiding the need to pre-merge onto a common grid. Takes precedence over eef_per_m in met if both are present. Assumed to adhere to pycontrails MetDataArray conventions.

  • dag (HorizontalDAG or None, default None) – Pre-built DAG. If None, a DAG is generated via Poisson-disk sampling along the great circle. The DAG origin and destination must agree with the airport coordinates.

  • cost_index (float, default 60.0) – Fuel-vs-time tradeoff in kg per minute. Higher values penalize time more, favoring faster (and more fuel-intensive) routes.

  • dollar_tonne_co2e (float, default 0.0) – Carbon price in US dollars per tonne (1000kg) of CO2-equivalent. A value of 0.0 disables the carbon cost term. If positive, either met must contain a eef_per_m variable or the eef parameter must be provided.

  • dollar_kg_fuel (float, default 1.0) – Fuel price in US dollars per kg. Only used to convert the carbon cost into the fuel-equivalent units of the objective function. Ignored if dollar_tonne_co2e is 0.0.

  • step_penalty_kg (float, default 0.0) – Cost in kg of fuel charged for changing flight level, on top of the maneuver’s own fuel and time, to discourage marginally beneficial steps. The initial climb and final descent are exempt. The penalty enters the objective only, not the aircraft mass, so reported fuel burn stays physical.

  • met_spacing_m (float, default 25_000.0) – Spacing in meters between met sample points along each edge.

  • flight_hours (int or None, default None) – Upper-bound flight duration in hours for met time window. If None, estimated from the aircraft type. Providing an explicit value decouples the met lookup from the aircraft, allowing the user to call the solve() method with a different aircraft type without re-initializing the optimizer.

  • allow_cooling_credit (bool, default False) – If True, negative EEF (cooling contrails) reduces cost when dollar_tonne_co2e is set. If False, negative EEF is clipped to zero in the cost function but still reported in the output flight. Only used if dollar_tonne_co2e is set.

  • avoidance_regions (list of polygon coordinate lists, or None) – Polygons to exclude from the search, defined as lists of (lon, lat) vertices. Edges intersecting any polygon are removed and the DAG is re-pruned.

  • **kwargs – Additional parameters for DAG generation if dag is None. Passed into HorizontalDAG.from_poisson.

classmethod from_flight(flight: Flight, *, met: MetDataset | Dataset | None = None, fl_profile: Dataset | None = None, aircraft_type: str | None = None, origin_icao: str | None = None, dest_icao: str | None = None, eef: DataArray | MetDataArray | None = None, altitude_ft: NDArray[floating] | None = None, cost_index: float = 60.0, dollar_tonne_co2e: float = 0.0, dollar_kg_fuel: float = 1.0, allow_cooling_credit: bool = False, use_flown_climb_descent: bool = False) Self[source]

Build a vertical-profile optimizer from a pycontrails.Flight trajectory.

The optimized flight follows the flight’s lateral path exactly, choosing flight level and Mach number along it via solve_track(). Weather comes from one of two sources:

  • met: raw gridded 4D met, used to interpolate the flight’s waypoints at each candidate flight level.

  • fl_profile: an already-interpolated (waypoint, altitude_ft) dataset carrying air_temperature, u_wind, v_wind, and optionally eef_per_m, aligned waypoint-for-waypoint with flight.

Parameters:
  • flight (pycontrails.Flight) – Trajectory supplying the lateral path, schedule, and (for use_flown_climb_descent) the flown altitude profile.

  • met (pycontrails.MetDataset or xarray.Dataset or None) – Gridded met to interpolate. Mutually exclusive with fl_profile.

  • fl_profile (xarray.Dataset or None) – Pre-interpolated per-waypoint met columns. Mutually exclusive with met.

  • aircraft_type (str or None) – PS model key. If None, taken from flight.attrs.

  • origin_icao (str or None) – ICAO codes. If None, taken from flight.attrs, else the nearest airport.

  • dest_icao (str or None) – ICAO codes. If None, taken from flight.attrs, else the nearest airport.

  • eef (xarray.DataArray or pycontrails.MetDataArray or None) – Effective energy forcing per meter, if supplied separately from met.

  • altitude_ft (numpy.ndarray or None) – Candidate flight levels in feet, used only with met. If None, the eastbound/westbound defaults from cruise_flight_levels() are used. With fl_profile the levels come from its altitude_ft coordinate.

  • cost_index (float, default 60.0) – Fuel-vs-time tradeoff in kg per minute.

  • dollar_tonne_co2e (float, default 0.0) – Carbon price per tonne CO2-equivalent. If positive, the met source must carry eef_per_m.

  • dollar_kg_fuel (float, default 1.0) – Fuel price per kg, converting carbon cost into fuel-equivalent units.

  • allow_cooling_credit (bool, default False) – If True, negative EEF reduces cost when dollar_tonne_co2e is set.

  • use_flown_climb_descent (bool, default False) – If True, the flown initial climb and final descent below the lowest candidate flight level are taken from flight unchanged, and only the cruise phase above it is optimized.

property eef_cost_factor: float

Compute the kg-fuel-equivalent cost per J of effective energy forcing.

property kind: str

Return the optimization variant.

Returns "track" for a fixed-path vertical-profile 2d optimizer (built via from_flight()), or "dag" for the full 4d lateral-plus-vertical DAG optimizer.

solve(n_iter: int = 3, cost_index: float | None = None, dollar_tonne_co2e: float | None = None, aircraft_type: str | None = None, payload: float | None = None, allow_cooling_credit: bool | None = None, step_penalty_kg: float | None = None) DAGResult[source]

Solve the trajectory optimization via shortest-path dynamic programming on the DAG.

This method iteratively re-solves the DAG to converge on takeoff mass.

  • Estimate trip fuel from great-circle distance and set initial takeoff mass as landing mass + estimated trip fuel, capped at MTOW

  • Solve the dynamic program to find the optimal path and trip fuel

  • Update takeoff mass as landing mass + trip fuel, capped at MTOW

  • Stop after the takeoff mass estimate converges or after n_iter iterations

Parameters:
  • n_iter (int, default 3) – Maximum number of mass-convergence iterations. Each iteration re-solves the full DP.

  • cost_index (float or None, default None) – If provided, updates self.cost_index before solving. This parameter is safe to vary between calls without rebuilding intermediate artifacts.

  • dollar_tonne_co2e (float or None, default None) – If provided, updates self.dollar_tonne_co2e before solving. Safe to vary between calls without rebuilding intermediate artifacts.

  • aircraft_type (str or None, default None) – If provided, updates self.aircraft_type, self.atyp, and self.mach_choices before solving. Safe to vary between calls without rebuilding the DAG or met lookup provided the met lookup was built with a sufficiently long flight_hours window to accommodate the new aircraft’s speed.

  • payload (float or None, default None) – Aircraft payload in kg if known. If None, this is estimated with pycontrails.

  • allow_cooling_credit (bool or None, default None) – If provided, updates self.allow_cooling_credit before solving.

Returns:

DP state, takeoff mass, trip fuel, payload, reserve fuel, and landing mass.

Return type:

DAGResult

reconstruct_path() tuple[NDArray[int64], NDArray[int64], NDArray[float32]][source]

Trace backpointers from destination to origin to recover the optimal path.

The returned arrays are ordered origin-first.

Geographic coordinates for each waypoint are available via self.dag.lon[path_h] and self.dag.lat[path_h]. Flight levels are self.fl_choices[path_fl_idx] for interior waypoints; the origin uses a sentinel index (len(fl_choices)) representing ground level.

Returns:

  • path_h (numpy.ndarray) – Horizontal node indices along the path. The first and last entries are origin and destination.

  • path_fl_idx (numpy.ndarray) – Flight level index at each node. The origin uses a sentinel ground_fl_idx = len(fl_choices); the destination uses the actual cruise FL from which the final descent begins.

  • path_mach (numpy.ndarray) – Cruise Mach number on each incoming leg. The first entry path_mach[0] is NaN (no incoming leg at the origin).

to_flight() Flight[source]

Return the optimal trajectory as a pycontrails.Flight.

When met data is available, waypoints are emitted at each edge sample point (~20 km spacing) with proper climb/descent altitude profiles and per-sample eef_per_m. Without met, falls back to one waypoint per DAG node.

The solve() method must be called first.

plot_met(altitude_ft: float | None = None, time: Timestamp | None = None, ax: TypeAliasForwardRef('cartopy.mpl.geoaxes.GeoAxes') | None = None, show_wind_quiver: bool = True, show_eef: bool = True, **kwargs) cartopy.mpl.geoaxes.GeoAxes[source]

Plot met data on DAG nodes for a given flight level and time.

Draws a wind quiver overlay. When eef_per_m is available in the met lookup, also draws a scatter plot colored by EEF.

Parameters:
  • altitude_ft (float or None) – Flight level in feet (e.g. 37000). Snaps to the nearest available level. If None, uses the first available level.

  • time (pandas.Timestamp or None) – Time to select. Snaps to the nearest available time step. If None, uses the first available time step.

  • ax (cartopy.mpl.geoaxes.GeoAxes or None) – Cartopy GeoAxes to plot on. If None, calls self.dag.plot() to create one.

  • **kwargs – Passed to ax.quiver.

Returns:

The axes with the met overlay.

Return type:

cartopy.mpl.geoaxes.GeoAxes

animate_solve(display_fl_idx: int | None = None, ax: TypeAliasForwardRef('cartopy.mpl.geoaxes.GeoAxes') | None = None) matplotlib.animation.FuncAnimation[source]

Re-run the DP with converged mass and return a wavefront animation.

solve() must be called first. This re-runs a single solve_dag pass with the converged amass_init, capturing wavefront snapshots. Nodes are colored by best_cost[:, display_fl_idx] for a single FL. Edges whose source node has been processed are shown in blue; remaining edges are shown muted. The optimal path is then revealed dest to origin, matching the backpointer reconstruction order.

Parameters:
  • display_fl_idx (int | None) – FL index into fl_choices to display costs for. If None, uses the FL the optimal path spends the most legs at.

  • ax (GeoAxes | None) – Cartopy GeoAxes to draw on. If None, a new figure is created.

Returns:

Wavefront animation with cost coloring and optimal path reveal.

Return type:

matplotlib.animation.FuncAnimation