4D Flight Optimization¶
Given an origin-destination pair, a takeoff time, and an aircraft type, we can search for the trajectory that minimizes an airline’s cost of flying it. We call this 4D flight optimization: the optimizer lays out a sequence of waypoints in latitude, longitude, altitude, and time, choosing a flight level and Mach number along the way to minimize a combined fuel + time + climate cost.
This notebook works through a complete end-to-end example of 4D flight optimization, for a single flight from Boston to Seattle.
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
from pycontrails import DiskCacheStore, Flight
from pycontrails.datalib.google_forecast import GoogleForecast
from pycontrails.physics import units
import contrailopt as co
Flight¶
Specify the flight to optimize: an origin-destination pair, a takeoff time, and an aircraft type. The cost index, which sets the fuel-versus-time tradeoff, is chosen later when solving.
origin = co.AirportCoords.from_icao("KBOS")
dest = co.AirportCoords.from_icao("KSEA")
aircraft_type = "B737"
takeoff_time = pd.Timestamp("2026-03-07T02:00:00")
Weather¶
We use ERA5 reanalysis from the public ARCO-ERA5 archive, which requires no credentials. It supplies the air temperature and winds that drive aircraft performance and ground speed during optimization.
Both the weather and the contrail forecast are easier to handle if we downselect in time first. This requires an upper bound on the flight duration, which we estimate below.
flight_hours = co.estimate_flight_hours(origin, dest)
t0 = takeoff_time
t1 = t0 + pd.Timedelta(hours=flight_hours)
data_vars = {
"temperature": "air_temperature",
"u_component_of_wind": "eastward_wind",
"v_component_of_wind": "northward_wind",
}
arco_store = "gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3"
arco_ds = xr.open_zarr(arco_store, chunks=None)
# Select required variables, rename for pycontrails, downselect in time
met = arco_ds[list(data_vars)].rename(data_vars).sel(time=slice(t0, t1))
Contrail forecast¶
The climate term needs gridded contrail forcing data analogous to weather: effective energy forcing, in joules per meter of flight. We use the Google Contrails Forecast API, which requires an API key.
gf = GoogleForecast(time=(t0, t1), variables="eeef_per_m", cachestore=DiskCacheStore())
eef = gf.open_metdataset()["expected_effective_energy_forcing"]
# Plot a slice at FL360 for demonstration purposes
eef_slice = eef.data.sel(level=units.ft_to_pl(36000), time=takeoff_time, method="nearest")
cmap = mpl.colors.LinearSegmentedColormap.from_list(
"RdBu_r_warm", mpl.colormaps["RdBu_r"](np.linspace(0.5, 1.0, 256))
)
_, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(10, 6))
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=1, edgecolor="gray")
eef_slice.plot.pcolormesh(
x="longitude",
y="latitude",
ax=ax,
cmap=cmap,
vmin=0,
vmax=1e8,
extend="max",
cbar_kwargs={"shrink": 0.6},
)
for ap in (origin, dest):
ax.scatter(*ap.coords, color="blue")
ax.text(
ap.longitude + 1.0,
ap.latitude - 2.0,
ap.icao_code,
color="blue",
)
ax.set_title(f"FL360 — {takeoff_time:%Y-%m-%d %H:%M} UTC")
ax.set_extent([-125, -66, 35, 55]);
Cost optimization¶
We build an Optimizer from the route, the weather, and the contrail forecast, then solve it. Constructing the optimizer is the most expensive step: it lays out the candidate waypoint graph and interpolates the gridded data onto its edges. Solving is comparatively cheap, and the cost settings can be varied between calls without rebuilding the full optimizer.
opt = co.Optimizer(
origin_icao=origin,
dest_icao=dest,
takeoff_time=takeoff_time,
met=met,
eef=eef,
aircraft_type=aircraft_type,
rng=np.random.default_rng(1234), # reproducibility
)
opt
Optimizer(KBOS -> KSEA, B737, 2026-03-07 02:00:00, 453 nodes, with met, unsolved)
# The opt.solve method returns low-level information involved
# during the algorithm runtime.
res = opt.solve(cost_index=30)
# Calling to_flight generates the actual flight plan as a pycontrails.Flight
fl_opt = opt.to_flight()
fl_opt.attrs["total_fuel_burn"] = res.trip_fuel # used later
fl_opt
| Attributes | |
|---|---|
| time | [2026-03-07 02:00:00, 2026-03-07 07:11:04.341796875] |
| longitude | [-122.30999755859375, -71.00800323486328] |
| latitude | [42.36199951171875, 48.39500045776367] |
| altitude | [6.099999904632568, 12192.0] |
| aircraft_type | B737 |
| total_fuel_burn | 12686.153343008664 |
| mach_number | air_temperature | u_wind | v_wind | node_index | sample_index | eef_per_m | longitude | latitude | time | altitude_ft | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.378369 | 219.641129 | 27.281603 | -15.864287 | 0 | 66 | 3.616092e+05 | -71.007896 | 42.361965 | 2026-03-07 02:00:00.000000000 | 20.000000 |
| 1 | 0.399647 | 219.698410 | 28.460867 | -15.031495 | -1 | 67 | 4.258720e+05 | -71.298264 | 42.389095 | 2026-03-07 02:01:56.291824341 | 3124.419922 |
| 2 | 0.422542 | 219.744339 | 29.093954 | -14.010304 | -1 | 68 | 6.315872e+05 | -71.588890 | 42.415489 | 2026-03-07 02:03:52.590499878 | 6229.022949 |
| 3 | 0.447196 | 219.777603 | 29.373415 | -12.620842 | -1 | 69 | 9.673769e+05 | -71.879753 | 42.441151 | 2026-03-07 02:05:48.883148193 | 9333.464844 |
| 4 | 0.652464 | 219.748032 | 29.501022 | -10.938803 | -1 | 70 | 1.301366e+06 | -72.170853 | 42.466076 | 2026-03-07 02:07:45.180206299 | 12438.024414 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 167 | 0.707188 | 208.125290 | 24.134476 | -37.823910 | -1 | 35054 | 2.791093e+05 | -121.041382 | 47.499542 | 2026-03-07 07:02:57.771484375 | 16860.978516 |
| 168 | 0.656193 | 208.136841 | 24.323151 | -36.512398 | -1 | 35055 | 6.980336e+05 | -121.358833 | 47.487957 | 2026-03-07 07:04:59.414062500 | 12753.935547 |
| 169 | 0.441585 | 208.155457 | 24.821608 | -35.423721 | -1 | 35056 | 1.327501e+06 | -121.676125 | 47.475487 | 2026-03-07 07:07:01.052734375 | 8647.125000 |
| 170 | 0.409878 | 208.132996 | 25.381109 | -34.412952 | -1 | 35057 | 2.194665e+06 | -121.993279 | 47.462151 | 2026-03-07 07:09:02.697265625 | 4540.074219 |
| 171 | 0.381111 | 208.081924 | 25.599571 | -33.295048 | 452 | 35058 | 2.486246e+06 | -122.310257 | 47.447937 | 2026-03-07 07:11:04.341796875 | 432.996094 |
172 rows × 11 columns
ax = opt.plot_met(altitude_ft=36000, show_eef=False)
fl_opt.plot(ax=ax, linewidth=2, color="blue")
ax.get_figure().set_size_inches(10, 6)
Sweep the cost index¶
The cost index sets what a minute of flight time is worth in kilograms of fuel. The sweep below is one example of a simulation you can run with a built optimizer. We re-solve at three values, reusing the optimizer from above. We compare the fuel vs time tradeoff and visualize resulting routes and vertical profiles.
The three routes diverge over a portion of the flight, and the step-climb schedules differ throughout the cruise.
fls = []
summaries = {}
for cost_index in [0, 20, 60]:
res = opt.solve(cost_index=cost_index)
fl = opt.to_flight()
fl.attrs["cost_index"] = cost_index
fls.append(fl)
summaries[f"cost index {cost_index}"] = {
"fuel_kg": res.trip_fuel,
"duration_minutes": fl.duration.total_seconds() / 60.0,
}
# Display a DataFrame for presentation only
pd.DataFrame.from_dict(summaries, orient="index").round(1)
| fuel_kg | duration_minutes | |
|---|---|---|
| cost index 0 | 12476.4 | 328.1 |
| cost index 20 | 12602.0 | 314.2 |
| cost index 60 | 12984.5 | 304.8 |
_, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=1, edgecolor="gray")
ax.add_feature(cfeature.STATES, linewidth=0.5, edgecolor="gray")
ax.gridlines(draw_labels=True, alpha=0.3)
for fl in fls:
ax.plot(
fl["longitude"],
fl["latitude"],
transform=ccrs.PlateCarree(),
linewidth=2,
label=f"cost index {fl.attrs['cost_index']}",
)
ax.set_extent([-125, -66, 35, 55])
ax.legend()
plt.show()
_, ax = plt.subplots(figsize=(10, 6))
for fl in fls:
fl.plot_profile(ax=ax, label=f"cost index {fl.attrs['cost_index']}")
ax.set_yticks(opt.fl_choices)
ax.set_ylim(28000, 42000)
ax.grid(True, axis="y")
ax.legend()
plt.show()
Contrail-aware optimization¶
Aviation warms the climate in two main ways: CO2 released by burning jet fuel, and radiative forcing from contrails. The two act on very different timescales. CO2 persists for centuries to millennia, altering the radiative balance continuously. Contrails last hours, but force far more strongly while they do. By historical coincidence, aviation’s warming to date is split roughly evenly between the two.
A climate-aware optimizer has to trade one against the other. We price contrail energy forcing as a CO2-equivalent mass and charge dollar_tonne_co2e for it, alongside the CO2 from the fuel the trajectory burns. The conversion rests on AGWP-style metrics. While imperfect, they make the tradeoff explicit.
res = opt.solve(cost_index=30, dollar_tonne_co2e=10)
fl_opt_con = opt.to_flight()
fl_opt_con.attrs["total_fuel_burn"] = res.trip_fuel # used later
fl_opt_con
| Attributes | |
|---|---|
| time | [2026-03-07 02:00:00, 2026-03-07 07:08:43.404296875] |
| longitude | [-122.30999755859375, -71.00800323486328] |
| latitude | [42.36199951171875, 48.39500045776367] |
| altitude | [6.099999904632568, 12192.0] |
| aircraft_type | B737 |
| total_fuel_burn | 12845.543968008657 |
| mach_number | air_temperature | u_wind | v_wind | node_index | sample_index | eef_per_m | longitude | latitude | time | altitude_ft | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.378369 | 225.109436 | 25.965668 | -15.457136 | 0 | 66 | 6.446373e+05 | -71.007896 | 42.361965 | 2026-03-07 02:00:00.000000000 | 20.000000 |
| 1 | 0.401191 | 225.056473 | 26.361698 | -14.018524 | -1 | 67 | 8.501729e+05 | -71.298264 | 42.389095 | 2026-03-07 02:01:56.077629089 | 3341.043945 |
| 2 | 0.425873 | 225.002960 | 26.555529 | -12.426508 | -1 | 68 | 1.154808e+06 | -71.588890 | 42.415489 | 2026-03-07 02:03:52.162078857 | 6662.283203 |
| 3 | 0.452594 | 224.980713 | 26.742624 | -10.656299 | -1 | 69 | 1.747330e+06 | -71.879753 | 42.441151 | 2026-03-07 02:05:48.240509033 | 9983.349609 |
| 4 | 0.662756 | 224.953217 | 27.087595 | -8.775844 | -1 | 70 | 2.283806e+06 | -72.170853 | 42.466076 | 2026-03-07 02:07:44.323364258 | 13304.542969 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 167 | 0.707188 | 208.128906 | 24.112701 | -37.967144 | -1 | 35054 | 2.535676e+05 | -121.041382 | 47.499542 | 2026-03-07 07:00:36.833984375 | 16860.978516 |
| 168 | 0.656193 | 208.137207 | 24.291971 | -36.644329 | -1 | 35055 | 6.258978e+05 | -121.358833 | 47.487957 | 2026-03-07 07:02:38.476562500 | 12753.935547 |
| 169 | 0.441585 | 208.153030 | 24.777712 | -35.540241 | -1 | 35056 | 1.175764e+06 | -121.676125 | 47.475487 | 2026-03-07 07:04:40.115234375 | 8647.125000 |
| 170 | 0.409878 | 208.128830 | 25.324558 | -34.512783 | -1 | 35057 | 1.940869e+06 | -121.993279 | 47.462151 | 2026-03-07 07:06:41.759765625 | 4540.074219 |
| 171 | 0.381111 | 208.077591 | 25.537754 | -33.381542 | 452 | 35058 | 2.217915e+06 | -122.310257 | 47.447937 | 2026-03-07 07:08:43.404296875 | 432.996094 |
172 rows × 11 columns
ax = fl_opt.plot_profile(label="Cost optimal")
fl_opt_con.plot_profile(ax=ax, label="Contrail aware")
eef_curtain = fl_opt_con.intersect_met_cross_section(eef)
eef_curtain.plot.pcolormesh(
x="time",
y="altitude_ft",
ax=ax,
cmap=cmap,
vmin=0,
vmax=1e8,
extend="max",
)
ax.set_yticks(opt.fl_choices)
ax.set_ylim(28000, 42000)
ax.grid(True, axis="y")
ax.legend()
plt.show()
ax = opt.plot_met(altitude_ft=36000, time=takeoff_time)
fl_opt.plot(ax=ax, linewidth=2, label="cost optimal")
fl_opt_con.plot(ax=ax, linewidth=2, label="contrail optimal")
ax.legend()
ax.get_figure().set_size_inches(10, 6);
The cost of contrail mitigation¶
The table below shows the tradeoff between estimated contrail warming, time, and fuel. Pricing contrails at $10 per tonne of CO2e removes the majority of the trajectory’s contrail forcing at the expense of a small fuel penalty. Because the contrail optimal route flies below a contrail forcing region rather than climbing through it, it flies in warmer air, and the total flight duration ends up slightly shorter than that of the cost-optimal route.
def summarize(fl: Flight) -> dict[str, float]:
"""Summarize the contrail forcing, fuel burn, and duration of a solved trajectory."""
seg_eef = fl["eef_per_m"] * fl.segment_length()
return {
"eef_TJ": np.nansum(seg_eef).item() / 1e12,
"fuel_kg": fl.attrs["total_fuel_burn"],
"duration_minutes": fl.duration.total_seconds() / 60.0,
}
# Display a DataFrame for presentation only
pd.DataFrame.from_dict(
{
"cost optimal": summarize(fl_opt),
"contrail aware": summarize(fl_opt_con),
},
orient="index",
).round(1)
| eef_TJ | fuel_kg | duration_minutes | |
|---|---|---|---|
| cost optimal | 51.6 | 12686.2 | 311.1 |
| contrail aware | 6.0 | 12845.5 | 308.7 |