Contrail-aware company route choice

This example uses the same six illustrative AMS-FCO company routes as the fuel-only route-choice example. It first establishes the minimum-fuel result, then repeats the optimization with a four-dimensional ERA5 persistent-contrail field and a blended fuel-plus-contrail objective.

The focus here is the climate-cost extension. See route_choice_optimization.ipynb for the detailed route abstractions, company waypoint-network search, and relationship with RAD.

1. Define the company-route catalogue

Each alternative is a RouteOption with origin, destination, and ordered intermediate waypoints. OpenTOP independently optimizes the speed and vertical profile for every route.

import warnings
from itertools import pairwise
from pathlib import Path

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import openap.casadi as oc
from openap.aero import kts

import numpy as np
import opentop as top
import pandas as pd
from pyproj import Geod

warnings.filterwarnings("ignore", message=".*wave drag.*")
plt.rcParams.update(
    {
        "figure.dpi": 120,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.grid": True,
        "grid.alpha": 0.22,
    }
)
COLORS = ("#0072B2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9")
GEOD = Geod(ellps="WGS84")

ROOT = Path.cwd()
if ROOT.name == "examples":
    ROOT = ROOT.parent
from opentop.plotting import (  # noqa: E402
    LINE_STYLES,
    MARKERS,
    apply_publication_style,
    plot_wind_vectors,
    style_map,
)

apply_publication_style()

AMS = (52.3086, 4.7639)
FCO = (41.8003, 12.2389)

routes = (
    top.RouteOption(
        "Company North",
        (AMS, (51.5, 7.2), (49.2, 8.0), (47.0, 10.8), (44.5, 11.0), FCO),
        metadata={"catalogue": "illustrative", "route_id": "NORTH-1"},
    ),
    top.RouteOption(
        "Company West",
        (AMS, (50.6, 3.7), (48.2, 5.5), (46.2, 7.5), (44.0, 9.8), FCO),
        metadata={"catalogue": "illustrative", "route_id": "WEST-1"},
    ),
    top.RouteOption(
        "Company East",
        (AMS, (51.0, 8.0), (49.2, 11.0), (47.0, 13.0), (44.5, 13.5), FCO),
        metadata={"catalogue": "illustrative", "route_id": "EAST-1"},
    ),
    top.RouteOption(
        "Company Central",
        (AMS, (50.7, 8.0), (48.5, 9.5), (46.2, 9.0), (43.8, 11.5), FCO),
        metadata={"catalogue": "illustrative", "route_id": "CENTRAL-1"},
    ),
    top.RouteOption(
        "Company Alpine",
        (AMS, (50.8, 3.8), (49.0, 4.5), (47.0, 6.0), (45.0, 7.5), (43.5, 10.2), FCO),
        metadata={"catalogue": "illustrative", "route_id": "ALPINE-1"},
    ),
    top.RouteOption(
        "Company Adriatic",
        (AMS, (51.1, 7.5), (49.5, 10.2), (47.5, 12.1), (45.5, 13.0), (43.5, 13.2), FCO),
        metadata={"catalogue": "illustrative", "route_id": "ADRIATIC-1"},
    ),
)

provided_routes = pd.DataFrame(
    {
        "route": [route.name for route in routes],
        "route_id": [route.metadata["route_id"] for route in routes],
        "waypoints": [len(route.waypoints) for route in routes],
        "distance_km": [route.distance_m / 1_000 for route in routes],
    }
).sort_values("distance_km")
provided_routes.round({"distance_km": 1})

	route	route_id	waypoints	distance_km
0	Company North	NORTH-1	6	1368.5
3	Company Central	CENTRAL-1	6	1375.1
1	Company West	WEST-1	6	1389.6
4	Company Alpine	ALPINE-1	7	1411.6
5	Company Adriatic	ADRIATIC-1	7	1417.3
2	Company East	EAST-1	6	1444.5

fig = plt.figure(figsize=(12, 4.8), constrained_layout=True)
ax_map = fig.add_subplot(1, 2, 1, projection=ccrs.LambertConformal(central_longitude=9))
ax_bar = fig.add_subplot(1, 2, 2)

ax_map.set_extent([2.0, 15.5, 40.0, 54.0], crs=ccrs.PlateCarree())
ax_map.add_feature(cfeature.LAND.with_scale("50m"), facecolor="#f3f1ec")
ax_map.add_feature(cfeature.OCEAN.with_scale("50m"), facecolor="#eaf2f8")
ax_map.add_feature(cfeature.BORDERS.with_scale("50m"), linewidth=0.45, edgecolor="0.55")
ax_map.coastlines(resolution="50m", linewidth=0.6, color="0.35")
for route, color, line_style, marker in zip(routes, COLORS, LINE_STYLES, MARKERS):
    latitudes, longitudes = zip(*route.waypoints)
    ax_map.plot(
        longitudes,
        latitudes,
        color=color,
        linestyle=line_style,
        marker=marker,
        linewidth=1.8,
        markersize=3.5,
        label=route.name,
        transform=ccrs.PlateCarree(),
    )
ax_map.scatter(
    [AMS[1], FCO[1]],
    [AMS[0], FCO[0]],
    s=45,
    marker="s",
    color="#222222",
    zorder=5,
    transform=ccrs.PlateCarree(),
)
ax_map.text(AMS[1] - 0.35, AMS[0] + 0.35, "AMS", transform=ccrs.PlateCarree())
ax_map.text(FCO[1] + 0.15, FCO[0] - 0.2, "FCO", transform=ccrs.PlateCarree())
ax_map.set_title("Provided route catalogue", loc="left", fontweight="bold")
ax_map.legend(loc="lower left", frameon=True, fontsize=7, ncol=2)

ordered = provided_routes.sort_values("distance_km", ascending=False)
color_by_name = dict(zip([route.name for route in routes], COLORS))
bars = ax_bar.barh(
    ordered.route,
    ordered.distance_km,
    color=[color_by_name[name] for name in ordered.route],
)
ax_bar.bar_label(bars, fmt="%.0f km", padding=4)
ax_bar.set_xlabel("Geodesic polyline distance (km)")
ax_bar.set_title(
    "Distance is useful for screening, not final selection",
    loc="left",
    fontweight="bold",
)
ax_bar.set_xlim(0, ordered.distance_km.max() * 1.18)
plt.show()

Reproducible ERA5 wind and contrail conditions

Both route stages use the same four-dimensional ERA5 weather sample. build_route_wind_sample.py contains the shared FastMeteo request, wind conversion, and persistent-contrail calculation. The cells below load checked-in caches for a fast, offline run; set REFRESH_FASTMETEO = True to regenerate both products through FastMeteo.

Persistent-contrail potential is calculated from the Schmidt-Appleman criterion and ice supersaturation, then spatially smoothed before construction of a CasADi bspline interpolant.

from runpy import run_path

WEATHER_HELPER = ROOT / "examples" / "data" / "build_route_wind_sample.py"
REFRESH_FASTMETEO = False

weather_helpers = run_path(str(WEATHER_HELPER))
wind = weather_helpers["load_fastmeteo_wind_sample"](refresh=REFRESH_FASTMETEO)
contrail_interpolant = weather_helpers["load_fastmeteo_contrail_sample"](
    refresh=REFRESH_FASTMETEO
)
print(f"Loaded {len(wind):,} FastMeteo/ERA5 wind samples")
print("Loaded 4-D persistent-contrail cost interpolant")

Loaded 2,688 FastMeteo/ERA5 wind samples
Loaded 4-D persistent-contrail cost interpolant

2. Establish the minimum-fuel baseline

The factory returns a fresh wind-enabled CompleteFlight for every option. Ordered route points constrain the horizontal path, while altitude, Mach, vertical rate, timing, and mass remain continuous decision variables. The baseline objective is fuel; no cost index is used.

def optimizer_factory():
    optimizer = top.CompleteFlight("A320", AMS, FCO, m0=0.85)
    optimizer.enable_wind(wind)
    return optimizer


optimization = top.RouteOptimizationConfig(
    objective="fuel",
    minimum_nodes=30,
    nodes_per_leg=2,
    waypoint_tolerance_m=15_000.0,
)
choice = top.optimize_routes(routes, optimizer_factory, config=optimization)

assert choice.best is not None
print(f"Selected: {choice.best.route.name}")
print(f"Optimized fuel: {choice.best.fuel_kg:,.0f} kg")

Selected: Company North
Optimized fuel: 5,009 kg

Wind field used by the optimization

Blue arrows show the eastward and northward ERA5 components at FL350. The minimum-fuel route is orange; the other company options are retained in gray for context. Arrow length is quantitative and referenced by the 20 m/s key.

snapshot = wind.query("h == 10668 and ts == 0")

fig = plt.figure(figsize=(9.0, 5.5), constrained_layout=True)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=9.0))
style_map(
    ax,
    extent=[2.0, 16.0, 40.0, 54.0],
    data_crs=ccrs.PlateCarree(),
    resolution="50m",
)

for route in routes:
    route_lats, route_lons = zip(*route.waypoints)
    color = "#D55E00" if route == choice.best.route else "#777777"
    width = 2.2 if route == choice.best.route else 0.9
    alpha = 1.0 if route == choice.best.route else 0.45
    ax.plot(
        route_lons,
        route_lats,
        color=color,
        linewidth=width,
        alpha=alpha,
        transform=ccrs.PlateCarree(),
        zorder=3,
    )

plot_wind_vectors(
    ax,
    snapshot,
    data_crs=ccrs.PlateCarree(),
    scale=320,
    width=0.0035,
    key_x=0.78,
)
ax.set_title(
    f"ERA5 wind at FL350 and optimized route: {choice.best.route.name} - "
    "1 May 2021, 08:00 UTC"
)

Text(0.5, 1.0, 'ERA5 wind at FL350 and optimized route: Company North - 1 May 2021, 08:00 UTC')

ranking = pd.DataFrame(
    {
        "route": [item.route.name for item in choice.optimized],
        "success": [item.success for item in choice.optimized],
        "fuel_kg": [item.fuel_kg for item in choice.optimized],
        "solve_s": choice.solve_seconds,
        "distance_km": [item.route.distance_m / 1_000 for item in choice.optimized],
        "status": [item.status for item in choice.optimized],
    }
).sort_values("fuel_kg")
ranking.round({"fuel_kg": 1, "solve_s": 1, "distance_km": 1})

	route	success	fuel_kg	solve_s	distance_km	status
0	Company North	True	5009.1	13.7	1368.5	Solve_Succeeded
3	Company Central	True	5025.4	12.0	1375.1	Solve_Succeeded
1	Company West	True	5087.6	14.2	1389.6	Solve_Succeeded
4	Company Alpine	True	5143.6	12.0	1411.6	Solve_Succeeded
5	Company Adriatic	True	5197.2	11.2	1417.3	Solve_Succeeded
2	Company East	True	5278.2	13.0	1444.5	Solve_Succeeded

The ranking above is based on fuel from the nonlinear aircraft-performance optimization. A graph or catalogue can rank by distance first to reduce the candidate count, but distance is not substituted for the final optimized metric.

def cumulative_distance_km(trajectory):
    distances = [0.0]
    points = list(zip(trajectory.latitude, trajectory.longitude))
    for (lat_a, lon_a), (lat_b, lon_b) in pairwise(points):
        distances.append(GEOD.inv(lon_a, lat_a, lon_b, lat_b)[2] / 1_000)
    return np.cumsum(distances)


fig = plt.figure(figsize=(13, 8.2), constrained_layout=True)
ax_map = fig.add_subplot(2, 2, 1, projection=ccrs.LambertConformal(central_longitude=9))
ax_alt = fig.add_subplot(2, 2, 2)
ax_mach = fig.add_subplot(2, 2, 3)
ax_fuel = fig.add_subplot(2, 2, 4)

ax_map.set_extent([2.0, 15.5, 40.0, 54.0], crs=ccrs.PlateCarree())
ax_map.add_feature(cfeature.LAND.with_scale("50m"), facecolor="#f3f1ec")
ax_map.add_feature(cfeature.OCEAN.with_scale("50m"), facecolor="#eaf2f8")
ax_map.add_feature(cfeature.BORDERS.with_scale("50m"), linewidth=0.45, edgecolor="0.55")
ax_map.coastlines(resolution="50m", linewidth=0.6, color="0.35")

for item, color, line_style in zip(choice.optimized, COLORS, LINE_STYLES):
    if not item.success:
        continue
    trajectory = item.trajectory
    distance = cumulative_distance_km(trajectory)
    selected = item is choice.best
    linewidth = 3.0 if selected else 1.7
    alpha = 1.0 if selected else 0.72
    label = item.route.name + (" (selected)" if selected else "")
    ax_map.plot(
        trajectory.longitude,
        trajectory.latitude,
        color=color,
        linestyle=line_style,
        linewidth=linewidth,
        alpha=alpha,
        label=label,
        transform=ccrs.PlateCarree(),
    )
    ax_alt.plot(
        distance,
        trajectory.altitude / 1_000,
        color=color,
        linestyle=line_style,
        linewidth=linewidth,
        alpha=alpha,
        label=label,
    )
    ax_mach.plot(
        distance,
        trajectory.mach,
        color=color,
        linestyle=line_style,
        linewidth=linewidth,
        alpha=alpha,
        label=label,
    )

ax_map.set_title("a  Optimized horizontal trajectories", loc="left", fontweight="bold")
ax_map.legend(loc="lower left", fontsize=7, frameon=True, ncol=2)
ax_alt.set(
    title="b  Optimized vertical profiles",
    xlabel="Distance along trajectory (km)",
    ylabel="Altitude (thousand ft)",
)
ax_mach.set(
    title="c  Optimized speed profiles",
    xlabel="Distance along trajectory (km)",
    ylabel="Mach",
)

fuel_order = ranking.sort_values("fuel_kg", ascending=False)
bars = ax_fuel.barh(
    fuel_order.route,
    fuel_order.fuel_kg,
    color=[color_by_name[name] for name in fuel_order.route],
)
ax_fuel.bar_label(bars, fmt="%.0f kg", padding=4)
ax_fuel.set(title="d  Final ranking by optimized fuel", xlabel="Fuel (kg)")
ax_fuel.set_xlim(0, fuel_order.fuel_kg.max() * 1.18)
plt.show()

3. Add a persistent-contrail formation cost

The grid value is a smoothed persistent-contrail potential between zero and one. Multiplying it by true airspeed and interval duration gives an expected persistent-contrail distance. We combine that distance with fuel using an explicit policy weight:

[ J = m_\mathrm{fuel} + w_\mathrm{contrail} L_\mathrm{persistent}, \qquad w_\mathrm{contrail}=20;\mathrm{kg_{eq},km^{-1}}. ]

The weight is intentionally visible and tunable: it is an operational trade-off parameter, not a universal conversion from contrail length to climate forcing. A production study should replace it with the operator's chosen climate metric or marginal cost.

objective_factory constructs the callable for each fresh optimizer because the grid cost uses that optimizer's map projection. ranking_metric="objective" then selects the lowest blended objective instead of silently falling back to the lowest-fuel route.

CONTRAIL_PENALTY_KG_PER_KM = 20.0


def fuel_contrail_objective_factory(optimizer):
    def objective(x, u, dt, **kwargs):
        true_airspeed = oc.aero.mach2tas(u[0], x[2])
        persistent_contrail_km = (
            optimizer.obj_grid_cost(
                x,
                u,
                dt,
                interpolant=kwargs["interpolant"],
                n_dim=4,
                time_dependent=True,
            )
            * true_airspeed
            * 1e-3
        )
        fuel_kg = optimizer.obj_fuel(x, u, dt)
        return fuel_kg + CONTRAIL_PENALTY_KG_PER_KM * persistent_contrail_km

    return objective


contrail_optimization = top.RouteOptimizationConfig(
    objective_factory=fuel_contrail_objective_factory,
    ranking_metric="objective",
    minimum_nodes=30,
    nodes_per_leg=2,
    waypoint_tolerance_m=15_000.0,
    trajectory_kwargs={
        "interpolant": contrail_interpolant,
        "n_dim": 4,
        "time_dependent": True,
        "auto_rescale_objective": True,
    },
)
contrail_choice = top.optimize_routes(
    routes,
    optimizer_factory,
    config=contrail_optimization,
)

assert contrail_choice.best is not None
print(f"Minimum-fuel selection: {choice.best.route.name}")
print(f"Contrail-aware selection: {contrail_choice.best.route.name}")

Minimum-fuel selection: Company North
Contrail-aware selection: Company Central

def persistent_contrail_km(trajectory):
    states = np.vstack(
        [
            trajectory.longitude,
            trajectory.latitude,
            trajectory.h,
            trajectory.ts,
        ]
    )
    potential = np.asarray(contrail_interpolant(states)).reshape(-1)
    segment_km = (
        trajectory.tas.to_numpy()[:-1] * kts * np.diff(trajectory.ts.to_numpy()) / 1_000
    )
    return float(np.sum(potential[:-1] * segment_km))


fuel_results = {item.route.name: item for item in choice.successful}
contrail_results = {item.route.name: item for item in contrail_choice.successful}
contrail_ranking = pd.DataFrame(
    {
        "route": [route.name for route in routes],
        "fuel_baseline_kg": [fuel_results[route.name].fuel_kg for route in routes],
        "fuel_contrail_aware_kg": [
            contrail_results[route.name].fuel_kg for route in routes
        ],
        "contrail_baseline_km": [
            persistent_contrail_km(fuel_results[route.name].trajectory)
            for route in routes
        ],
        "contrail_aware_km": [
            persistent_contrail_km(contrail_results[route.name].trajectory)
            for route in routes
        ],
        "blended_objective": [
            contrail_results[route.name].objective_value for route in routes
        ],
    }
).sort_values("blended_objective")
contrail_ranking.round(1)

	route	fuel_baseline_kg	fuel_contrail_aware_kg	contrail_baseline_km	contrail_aware_km	blended_objective
3	Company Central	5025.4	5128.5	88.1	81.3	6758.1
1	Company West	5087.6	5150.3	89.5	85.5	6867.5
0	Company North	5009.1	5157.1	89.5	86.4	6882.0
4	Company Alpine	5143.6	5221.0	90.3	87.9	6980.0
5	Company Adriatic	5197.2	5655.1	104.2	70.2	7064.4
2	Company East	5278.2	5672.4	105.3	79.3	7265.3

The map below shows a horizontal slice of the same four-dimensional field used by the optimizer. The trajectory optimization evaluates the full field at every collocation point, so it can avoid persistent-contrail regions vertically as well as by selecting a different company route.

grid_longitude = np.linspace(2.0, 16.0, 100)
grid_latitude = np.linspace(40.0, 54.0, 100)
grid_lon, grid_lat = np.meshgrid(grid_longitude, grid_latitude)
field_input = np.vstack(
    [
        grid_lon.ravel(),
        grid_lat.ravel(),
        np.full(grid_lon.size, 35_000 * 0.3048),
        np.full(grid_lon.size, 3_600.0),
    ]
)
contrail_field = np.asarray(contrail_interpolant(field_input)).reshape(grid_lon.shape)

fig = plt.figure(figsize=(13, 8.3), constrained_layout=True)
ax_map = fig.add_subplot(2, 2, 1, projection=ccrs.LambertConformal(central_longitude=9))
ax_altitude = fig.add_subplot(2, 2, 2)
ax_fuel = fig.add_subplot(2, 2, 3)
ax_contrail = fig.add_subplot(2, 2, 4)
style_map(
    ax_map,
    extent=[2.0, 16.0, 40.0, 54.0],
    data_crs=ccrs.PlateCarree(),
    resolution="50m",
)

field = ax_map.contourf(
    grid_lon,
    grid_lat,
    np.clip(contrail_field, 0.0, 1.0),
    levels=np.linspace(0.0, 0.25, 11),
    cmap="viridis",
    alpha=0.72,
    transform=ccrs.PlateCarree(),
)
fig.colorbar(
    field,
    ax=ax_map,
    orientation="horizontal",
    pad=0.03,
    label="Persistent-contrail potential at FL350, +1 h",
)

for route in routes:
    route_latitude, route_longitude = zip(*route.waypoints)
    ax_map.plot(
        route_longitude,
        route_latitude,
        color="#777777",
        linewidth=0.8,
        alpha=0.45,
        transform=ccrs.PlateCarree(),
    )

selected_cases = (
    ("Minimum fuel", choice.best, COLORS[0], "-"),
    ("Fuel + contrail", contrail_choice.best, COLORS[1], "--"),
)
for label, result, color, line_style in selected_cases:
    trajectory = result.trajectory
    distance = cumulative_distance_km(trajectory)
    ax_map.plot(
        trajectory.longitude,
        trajectory.latitude,
        color=color,
        linestyle=line_style,
        linewidth=2.8,
        label=f"{label}: {result.route.name}",
        transform=ccrs.PlateCarree(),
    )
    ax_altitude.plot(
        distance,
        trajectory.altitude / 1_000,
        color=color,
        linestyle=line_style,
        linewidth=2.4,
        label=label,
    )

ax_map.set_title("a  Contrail field and selected trajectories", loc="left")
ax_map.legend(loc="lower left", fontsize=7)
ax_altitude.set(
    title="b  Selected vertical profiles",
    xlabel="Distance along trajectory (km)",
    ylabel="Altitude (thousand ft)",
)
ax_altitude.legend()

plot_order = contrail_ranking.sort_values(
    "blended_objective", ascending=True
).reset_index(drop=True)
y = np.arange(len(plot_order))
bar_height = 0.36
ax_fuel.barh(
    y - bar_height / 2,
    plot_order.fuel_baseline_kg,
    height=bar_height,
    color=COLORS[0],
    label="Minimum-fuel solve",
)
ax_fuel.barh(
    y + bar_height / 2,
    plot_order.fuel_contrail_aware_kg,
    height=bar_height,
    color=COLORS[1],
    label="Fuel + contrail solve",
)
ax_fuel.set(
    title="c  Fuel trade-off by route",
    xlabel="Fuel (kg)",
    yticks=y,
    yticklabels=plot_order.route,
)
ax_fuel.legend(loc="upper center", bbox_to_anchor=(0.5, -0.14), ncol=2)

ax_contrail.barh(
    y - bar_height / 2,
    plot_order.contrail_baseline_km,
    height=bar_height,
    color=COLORS[0],
    label="Minimum-fuel solve",
)
ax_contrail.barh(
    y + bar_height / 2,
    plot_order.contrail_aware_km,
    height=bar_height,
    color=COLORS[1],
    label="Fuel + contrail solve",
)
ax_contrail.set(
    title="d  Expected persistent-contrail distance",
    xlabel="Persistent-contrail distance (km)",
    yticks=y,
    yticklabels=plot_order.route,
)
plt.show()

Takeaways

• The fuel-only and contrail-aware cases use the same routes, wind, aircraft, and optimization API.
• objective_factory creates the projection-aware grid objective for each route solve.
• ranking_metric="objective" selects by the blended cost rather than fuel alone.
• The contrail penalty is an explicit policy parameter and should be replaced with the climate metric appropriate to the study.