Multiple-aircraft separation optimization

This notebook studies crossing scenarios with 2, 3, 4, 6, and 8 cruise flights. It is divided into two sections:

1. Altitude free: the original scenarios, where aircraft may use both lateral and vertical avoidance.
2. Common fixed altitude: level aircraft constrained to the same optimized altitude, so conflicts must be avoided horizontally.

Each scenario uses one CasADi Opti problem and the default high-order separation approximation:

(dx / 5 NM)^2 + (dy / 5 NM)^2 + (dh / 1,000 ft)^8 >= 1.3.

from dataclasses import asdict

import matplotlib.pyplot as plt
import openap

import numpy as np
import opentop as top
import pandas as pd

%matplotlib inline

Shared solver helper

The helper below solves one family of scenarios and returns a compact five-row summary. Detailed pairwise reports remain available separately instead of expanding the displayed table.

scenario_sizes = (2, 3, 4, 6, 8)


def solve_scenarios(scenarios, *, common_altitude=False):
    results = {}
    reports_by_count = {}
    rows = []
    for count, flights in scenarios.items():
        result = top.MultiAircraft(
            flights,
            max_iter=2_000,
            common_altitude=common_altitude,
        ).trajectory()
        results[count] = result
        reports = pd.DataFrame(asdict(report) for report in result.pair_reports)
        reports_by_count[count] = reports
        worst = min(result.pair_reports, key=lambda report: report.minimum_metric)
        row = {
            "aircraft": count,
            "passed": result.success and result.separation_success,
            "worst pair": f"{worst.first} / {worst.second}",
            "min metric": worst.minimum_metric,
            "solve (s)": result.solve_time_s,
        }
        if common_altitude:
            first_trajectory = next(iter(result.trajectories.values()))
            row["shared altitude (ft)"] = first_trajectory.altitude.iloc[0]
        rows.append(row)

    summary = pd.DataFrame(rows).set_index("aircraft")
    return (
        results,
        reports_by_count,
        summary.round(
            {
                "min metric": 3,
                "solve (s)": 1,
                "shared altitude (ft)": 0,
            }
        ),
    )

1. Altitude-free scenarios

This section preserves the earlier notebook setup. Routes use bearings around the full 360 degrees, so some aircraft fly reciprocal tracks. The 2-, 3-, and 4-aircraft cases start together. The 6- and 8-aircraft cases use two arrival waves separated by three minutes. Altitude is free, allowing the eighth-order vertical term to contribute to conflict avoidance.

def build_free_altitude_scenario(count, *, nodes=8):
    center_lat, center_lon = 51.0, 7.0
    half_route_m = 150_000.0
    wave_size = count if count <= 4 else count // 2
    flights = []
    for index in range(count):
        bearing = 360.0 * index / count
        origin = openap.aero.latlon(center_lat, center_lon, half_route_m, bearing)
        destination = openap.aero.latlon(
            center_lat, center_lon, half_route_m, (bearing + 180.0) % 360.0
        )
        optimizer = top.Cruise(
            "A320",
            (float(origin[0]), float(origin[1])),
            (float(destination[0]), float(destination[1])),
            m0=0.85,
        )
        optimizer.setup(nodes=nodes, max_iter=2_000)
        flights.append(
            top.FlightSpec(
                f"AC{index + 1}",
                optimizer,
                start_time=180.0 * (index // wave_size),
            )
        )
    return flights


free_scenarios = {
    count: build_free_altitude_scenario(count) for count in scenario_sizes
}
{count: len(flights) for count, flights in free_scenarios.items()}

{2: 2, 3: 3, 4: 4, 6: 6, 8: 8}

Altitude-free results

The minimum metric is densely verified after each solve; the required value is 1.3.

free_results, free_pair_reports, free_summary = solve_scenarios(free_scenarios)
free_summary

	passed	worst pair	min metric	solve (s)
aircraft
2	True	AC1 / AC2	1.324	14.6
3	True	AC1 / AC3	1.311	17.9
4	True	AC2 / AC3	1.699	87.7
6	True	AC5 / AC6	3.133	113.2
8	True	AC7 / AC8	1.307	122.0

assert all(result.success for result in free_results.values())
assert all(result.separation_success for result in free_results.values())

Altitude-free plots

Each scenario has a ground-track and altitude panel. Reciprocal tracks and similar altitude profiles may overlap exactly, so small perpendicular-track and 100 ft altitude display offsets are added. These offsets are visualization-only; optimization and separation verification use the original trajectories.

def plot_free_altitude_scenario(count, result):
    fig, (route_ax, altitude_ax) = plt.subplots(1, 2, figsize=(12, 4.5))
    trajectory_count = len(result.trajectories)
    for index, (flight_id, trajectory) in enumerate(result.trajectories.items()):
        lane = index - (trajectory_count - 1) / 2
        d_lon = float(trajectory.longitude.iloc[-1] - trajectory.longitude.iloc[0])
        d_lat = float(trajectory.latitude.iloc[-1] - trajectory.latitude.iloc[0])
        route_norm = max(np.hypot(d_lon, d_lat), 1e-12)
        route_offset_deg = 0.012 * lane
        display_longitude = trajectory.longitude - d_lat / route_norm * route_offset_deg
        display_latitude = trajectory.latitude + d_lon / route_norm * route_offset_deg
        display_altitude = trajectory.altitude + 100.0 * lane

        route_ax.plot(display_longitude, display_latitude, marker="o", label=flight_id)
        altitude_ax.plot(
            trajectory.absolute_ts / 60.0,
            display_altitude,
            marker="o",
            label=flight_id,
        )

    route_ax.set(
        title="Ground tracks (display offsets)",
        xlabel="longitude (deg)",
        ylabel="latitude (deg)",
    )
    route_ax.set_aspect("equal", adjustable="datalim")
    altitude_ax.set(
        title="Altitude profiles (+100 ft display lanes)",
        xlabel="absolute time (min)",
        ylabel="display altitude (ft)",
    )
    for axis in (route_ax, altitude_ax):
        axis.grid(ls=":", alpha=0.7)
        axis.legend()
    fig.suptitle(f"{count} aircraft — altitude free")
    fig.tight_layout()
    plt.show()


for count, result in free_results.items():
    plot_free_altitude_scenario(count, result)

2. Common fixed-altitude scenarios

This separate horizontal-avoidance study uses the same aircraft counts. Every child optimizer calls fix_cruise_altitude(), making each trajectory level, and common_altitude=True equates those levels. The shared flight level remains an optimization variable rather than being prescribed.

Unlike the original altitude-free set, these scenarios use non-duplicated route lines and two-aircraft arrival waves separated by three minutes. Aircraft within each wave still encounter one another simultaneously. This avoids an unnecessarily pathological all-aircraft-at-one-point problem while preserving genuine horizontal coordination.

def build_common_altitude_scenario(count, *, nodes=8):
    center_lat, center_lon = 51.0, 7.0
    half_route_m = 150_000.0
    flights = []
    for index in range(count):
        bearing = 180.0 * index / count
        origin = openap.aero.latlon(center_lat, center_lon, half_route_m, bearing)
        destination = openap.aero.latlon(
            center_lat, center_lon, half_route_m, (bearing + 180.0) % 360.0
        )
        optimizer = top.Cruise(
            "A320",
            (float(origin[0]), float(origin[1])),
            (float(destination[0]), float(destination[1])),
            m0=0.85,
        )
        optimizer.setup(nodes=nodes, max_iter=2_000)
        optimizer.fix_cruise_altitude()
        flights.append(
            top.FlightSpec(
                f"AC{index + 1}",
                optimizer,
                start_time=180.0 * (index // 2),
            )
        )
    return flights


common_altitude_scenarios = {
    count: build_common_altitude_scenario(count) for count in scenario_sizes
}
{count: len(flights) for count, flights in common_altitude_scenarios.items()}

{2: 2, 3: 3, 4: 4, 6: 6, 8: 8}

Common-altitude results

The table includes the optimized shared altitude. Within each scenario, every aircraft has this same constant altitude and therefore zero vertical separation.

common_results, common_pair_reports, common_summary = solve_scenarios(
    common_altitude_scenarios, common_altitude=True
)
common_summary

	passed	worst pair	min metric	solve (s)	shared altitude (ft)
aircraft
2	True	AC1 / AC2	1.310	16.1	33762.0
3	True	AC1 / AC2	1.308	19.6	33754.0
4	True	AC3 / AC4	1.306	47.2	33708.0
6	True	AC1 / AC2	1.308	36.6	33649.0
8	True	AC7 / AC8	1.306	85.8	33581.0

assert all(result.success for result in common_results.values())
assert all(result.separation_success for result in common_results.values())
for result in common_results.values():
    altitudes = [
        trajectory.altitude.to_numpy() for trajectory in result.trajectories.values()
    ]
    assert np.ptp(np.concatenate(altitudes)) <= 1e-6

Horizontal conflict-avoidance plots

These figures show actual ground tracks without display offsets. The second panel includes every aircraft pair with an overlapping time window. Non-worst pairs are thin gray curves, while the densely verified worst pair is highlighted. With zero vertical difference, the high-order constraint reduces to a horizontal threshold of 5 NM * sqrt(1.3), approximately 5.7 NM.

def pair_horizontal_profile(result, report, samples=300):
    first = result.trajectories[report.first]
    second = result.trajectories[report.second]
    start = max(first.absolute_ts.iloc[0], second.absolute_ts.iloc[0])
    end = min(first.absolute_ts.iloc[-1], second.absolute_ts.iloc[-1])
    if end <= start:
        return np.array([]), np.array([])
    times = np.linspace(start, end, samples)
    dx = np.interp(times, first.absolute_ts, first.x) - np.interp(
        times, second.absolute_ts, second.x
    )
    dy = np.interp(times, first.absolute_ts, first.y) - np.interp(
        times, second.absolute_ts, second.y
    )
    return times, np.hypot(dx, dy) / 1852.0


def plot_common_altitude_scenario(count, result):
    fig, (route_ax, distance_ax) = plt.subplots(1, 2, figsize=(12, 4.5))
    for flight_id, trajectory in result.trajectories.items():
        route_ax.plot(
            trajectory.longitude,
            trajectory.latitude,
            marker="o",
            label=flight_id,
        )

    worst_report = min(result.pair_reports, key=lambda report: report.minimum_metric)
    other_label_available = True
    for report in result.pair_reports:
        if report == worst_report or not report.overlaps:
            continue
        times, horizontal_nm = pair_horizontal_profile(result, report)
        if not len(times):
            continue
        distance_ax.plot(
            times / 60.0,
            horizontal_nm,
            color="0.6",
            alpha=0.55,
            lw=1.0,
            label="other pairs" if other_label_available else "_nolegend_",
        )
        other_label_available = False

    times, horizontal_nm = pair_horizontal_profile(result, worst_report)
    config = top.SeparationConfig()
    required_nm = config.horizontal_m * np.sqrt(config.minimum_metric) / 1852.0
    distance_ax.plot(
        times / 60.0,
        horizontal_nm,
        color="tab:blue",
        lw=2.2,
        label=f"worst: {worst_report.first} / {worst_report.second}",
    )
    distance_ax.axhline(required_nm, color="tab:red", ls="--", label="required")
    distance_ax.scatter(
        worst_report.worst_time / 60.0,
        worst_report.horizontal_m / 1852.0,
        color="black",
        zorder=3,
        label="verified minimum",
    )

    route_ax.set(
        title="Ground tracks (actual)",
        xlabel="longitude (deg)",
        ylabel="latitude (deg)",
    )
    route_ax.set_aspect("equal", adjustable="datalim")
    distance_ax.set(
        title="All-pair horizontal separation",
        xlabel="absolute time (min)",
        ylabel="horizontal distance (NM)",
    )
    for axis in (route_ax, distance_ax):
        axis.grid(ls=":", alpha=0.7)
        axis.legend()
    shared_altitude = next(iter(result.trajectories.values())).altitude.iloc[0]
    fig.suptitle(f"{count} aircraft — common altitude {shared_altitude:,.0f} ft")
    fig.tight_layout()
    plt.show()


for count, result in common_results.items():
    plot_common_altitude_scenario(count, result)

Detailed pair diagnostics

Use pair_details() only when pair-level values are needed. This keeps the main result tables compact.

def pair_details(altitude_mode, count):
    report_sets = {
        "free": free_pair_reports,
        "common": common_pair_reports,
    }
    columns = [
        "first",
        "second",
        "minimum_metric",
        "required_metric",
        "margin",
        "worst_time",
        "horizontal_m",
        "vertical_m",
    ]
    return report_sets[altitude_mode][count][columns].round(3)


print('Examples: pair_details("free", 8), pair_details("common", 8)')

Examples: pair_details("free", 8), pair_details("common", 8)

Optional separation customization

Most users can rely on SeparationConfig() defaults. To study a different standard, pass a customized configuration to MultiAircraft; keep an even vertical_power to preserve a symmetric conflict region.

custom_separation = top.SeparationConfig(
    horizontal_m=6 * 1852.0,
    vertical_m=1_000 * 0.3048,
    vertical_power=8,
)

# Example:
# custom_result = top.MultiAircraft(
#     free_scenarios[2],
#     separation=custom_separation,
#     max_iter=2_000,
# ).trajectory()
custom_separation

SeparationConfig(horizontal_m=11112.0, vertical_m=304.8, vertical_power=8, minimum_metric=1.3, watch_factor=1.5, detect_dt_s=15.0, constraint_dt_s=10.0, encounter_buffer_s=35.0, verification_dt_s=5.0, max_refinements=2, verification_tolerance=0.001, constraint_buffer=0.01)