Six-aircraft crossing conflict

Six identical A320s start at the same time on approximately 300 km routes crossing near 51 N, 7 E. Every flight enters and exits at 30,000 ft. Independent fuel-optimal trajectories create a simultaneous conflict involving all 15 aircraft pairs. Two coordinated cases use the same geometry and boundary conditions:

1. Free interior altitude: aircraft may climb and descend within 27,000-33,000 ft, return to the fixed exit altitude, and adjust their horizontal tracks and flight durations.
2. Fixed altitude throughout: aircraft remain at 30,000 ft and resolve the conflict through horizontal-track and timing adjustments.

Each coordinated case uses one shared CasADi Opti problem and minimizes total fuel burn. Both use the default separation metric and threshold:

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

scenario_sizes = (6,) runs only the six-aircraft case shown here. Other counts, such as 2, 3, 4, or 8, can be added for additional solves; their results are not included in this execution.

from dataclasses import asdict
from itertools import combinations

import matplotlib.pyplot as plt
import openap

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

%matplotlib inline
from opentop.plotting import (
    OKABE_ITO,
    apply_publication_style,
)

apply_publication_style()

Shared solver settings and summary

Both cases use identical aircraft masses (m0=0.80), simultaneous entry at zero seconds, and fixed 30,000 ft endpoint altitudes. Each aircraft has ten collocation intervals. The joint solver uses exact Hessians and permits up to 2,000 iterations per solve.

Separation constraints are initially sampled every 30 seconds. The default verification checks the collocation trajectories every five seconds and permits up to two refinement rounds. The summary reports combined solver/separation success, the pair with the smallest verified metric, that metric, and solver time. The fixed-altitude case also reports its prescribed flight level. Full pair reports are retained separately.

scenario_sizes = (6,)
ENDPOINT_ALTITUDE_FT = 30_000.0
CONFLICT_CENTER = (51.0, 7.0)


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,
            separation=top.SeparationConfig(constraint_dt_s=30.0),
            ipopt_kwargs={"hessian_approximation": "exact"},
        ).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. Fixed endpoints, free interior altitude

The six origin bearings are 0, 30, 60, 90, 120, and 150 degrees from the shared center. Origins and destinations lie 150 km from the center on opposite bearings, giving distinct crossing tracks. All flights start at zero seconds.

EndpointCruise fixes the initial and final altitude bounds to 30,000 ft. Interior altitude remains bounded between 27,000 and 33,000 ft. Calling allow_cruise_descent() enables descent as well as climb; the aircraft need not remain level. Horizontal motion uses the route's projected bounding box expanded by 30 km in each coordinate direction. The center is not imposed as a waypoint, so coordinated tracks can avoid it.

The baseline check first solves every flight independently. It linearly interpolates the returned trajectories onto 2,000 shared time samples and requires all 15 pair metrics to fall below 1.3 at the same sample. The displayed conflict time is the first sampled instant satisfying that condition, not necessarily the instant of closest approach.

class EndpointCruise(top.Cruise):
    """Fix entry and exit altitudes while leaving interior altitude free."""

    def init_conditions(self, **kwargs):
        super().init_conditions(**kwargs)
        endpoint_h = ENDPOINT_ALTITUDE_FT * openap.aero.ft
        self.x_0_lb[2] = self.x_0_ub[2] = endpoint_h
        self.x_f_lb[2] = self.x_f_ub[2] = endpoint_h


def build_free_altitude_scenario(count, *, nodes=10):
    flights = []
    for index in range(count):
        bearing = 180.0 * index / count
        origin = openap.aero.latlon(*CONFLICT_CENTER, 150_000.0, bearing)
        destination = openap.aero.latlon(
            *CONFLICT_CENTER, 150_000.0, (bearing + 180.0) % 360.0
        )
        optimizer = EndpointCruise(
            "A320",
            (float(origin[0]), float(origin[1])),
            (float(destination[0]), float(destination[1])),
            m0=0.80,
            h_min=27_000.0 * openap.aero.ft,
            h_max=33_000.0 * openap.aero.ft,
        )
        optimizer.allow_cruise_descent()
        optimizer.setup(nodes=nodes, max_iter=2_000)
        flights.append(
            top.FlightSpec(
                f"AC{index + 1}",
                optimizer,
                start_time=0.0,
                options={"route_margin_m": 30_000.0},
            )
        )
    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()}

{6: 6}

def independent_conflict_check(scenarios):
    projection = Proj(
        proj="aeqd", lat_0=CONFLICT_CENTER[0], lon_0=CONFLICT_CENTER[1], datum="WGS84"
    )
    config = top.SeparationConfig()
    rows = []
    for count, flights in scenarios.items():
        independent = [
            flight.optimizer.trajectory(**dict(flight.options)) for flight in flights
        ]
        times = np.linspace(0.0, min(df.ts.iloc[-1] for df in independent), 2_000)
        positions = []
        for df in independent:
            x, y = projection(
                np.interp(times, df.ts, df.longitude),
                np.interp(times, df.ts, df.latitude),
            )
            positions.append((x, y, np.interp(times, df.ts, df.h)))
        metrics = np.array(
            [
                ((a[0] - b[0]) / config.horizontal_m) ** 2
                + ((a[1] - b[1]) / config.horizontal_m) ** 2
                + ((a[2] - b[2]) / config.vertical_m) ** config.vertical_power
                for a, b in combinations(positions, 2)
            ]
        )
        simultaneous = np.all(metrics < config.minimum_metric, axis=0)
        assert np.any(simultaneous), "Every pair must share one baseline conflict"
        rows.append(
            {
                "aircraft": count,
                "conflicting pairs": len(metrics),
                "shared conflict time (s)": times[np.flatnonzero(simultaneous)[0]],
            }
        )
    return pd.DataFrame(rows).set_index("aircraft").round(1)


independent_conflict_check(free_scenarios)

	conflicting pairs	shared conflict time (s)
aircraft
6	15	612.0

Free-interior-altitude results

The joint solve coordinates all six aircraft. A passing result requires solver success and separation verification within the configured tolerance. The assertions also check the aircraft/pair counts and confirm that every returned entry and exit altitude is 30,000 ft within 0.001 ft.

free_results, free_pair_reports, free_summary = solve_scenarios(free_scenarios)
free_summary

	passed	worst pair	min metric	solve (s)
aircraft
6	True	AC3 / AC5	1.307	193.8

def check_endpoints(results):
    for count, result in results.items():
        assert result.success and result.separation_success
        assert len(result.trajectories) == count
        assert len(result.pair_reports) == count * (count - 1) // 2
        for trajectory in result.trajectories.values():
            np.testing.assert_allclose(
                trajectory.altitude.iloc[[0, -1]],
                ENDPOINT_ALTITUDE_FT,
                atol=1e-3,
                rtol=0,
            )


check_endpoints(free_results)

Ground tracks and altitude profiles

The panels show the coordinated tracks and altitude profiles at returned trajectory nodes, with no display offsets. The dotted line marks the prescribed 30,000 ft entry and exit altitude. The plot label “altitude free” refers only to the interior profile within its altitude bounds; endpoint altitudes remain fixed. Flight durations may differ even though all aircraft start together.

def plot_free_altitude_scenario(count, result):
    fig, (route_ax, altitude_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
        )
        altitude_ax.plot(
            trajectory.absolute_ts / 60.0,
            trajectory.altitude,
            marker="o",
            label=flight_id,
        )
    altitude_ax.axhline(ENDPOINT_ALTITUDE_FT, color="0.5", ls=":")
    route_ax.set(
        title="Ground tracks",
        xlabel="longitude (deg)",
        ylabel="latitude (deg)",
    )
    route_ax.set_aspect("equal", adjustable="datalim")
    altitude_ax.set(
        title="Altitude profiles",
        xlabel="absolute time (min)",
        ylabel="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. Fixed altitude throughout

This comparison reuses the same scenario builder, including routes, aircraft masses, simultaneous entry times, and fixed endpoint altitudes. Each optimizer additionally calls fix_cruise_altitude(), which holds the whole trajectory level at 30,000 ft. This level constraint prevents climbs and descents even though the builder enables cruise descent. common_altitude=True explicitly equates the aircraft flight levels, which are already fixed by their endpoints.

The independent baseline is checked again for a simultaneous all-pair conflict. The coordinated solve can adjust horizontal tracks and flight durations; vertical avoidance is unavailable.

def build_common_altitude_scenario(count, *, nodes=10):
    flights = build_free_altitude_scenario(count, nodes=nodes)
    for flight in flights:
        flight.optimizer.fix_cruise_altitude()
    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()}

{6: 6}

Fixed-altitude baseline and coordinated results

The first table checks the independent baseline, using the same sampled conflict test as the first case. The second summarizes the coordinated solve. Its shared altitude is prescribed at 30,000 ft, not optimized freely. Assertions check solve/separation success, endpoint altitudes, and a common constant altitude across all returned trajectories.

independent_conflict_check(common_altitude_scenarios)

	conflicting pairs	shared conflict time (s)
aircraft
6	15	618.4

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
6	True	AC3 / AC4	1.307	104.7	30000.0

check_endpoints(common_results)
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 separation plots

The left panel shows actual coordinated ground tracks. The right panel shows all 15 pairwise horizontal-distance curves, calculated by linear interpolation between returned nodes at 300 time samples per pair. Gray curves show the other pairs; the highlighted curve corresponds to the worst pair identified by the solver's dense verification.

The black marker shows that pair's verified minimum, evaluated using the collocation trajectory. It can differ slightly from the displayed interpolated curve. With zero vertical separation, the required metric of 1.3 corresponds to a horizontal distance 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=OKABE_ITO[0],
        lw=2.2,
        label=f"worst: {worst_report.first} / {worst_report.second}",
    )
    distance_ax.axhline(required_nm, color=OKABE_ITO[1], 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("free", 6) for the free-interior-altitude case or pair_details("common", 6) for the fixed-altitude case. Each returns all 15 pair reports, including the minimum metric, required metric, margin, worst time in seconds, and horizontal and vertical distances in meters at that time. Other counts are available only if they were included in scenario_sizes and solved.

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", 6), pair_details("common", 6)')

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

Optional separation customization

The following cell only constructs and displays an alternative configuration; it does not rerun either case above. It changes the horizontal scale from 5 NM to 6 NM, retains the 1,000 ft vertical scale and eighth-power term, and uses defaults for the remaining settings, including the metric threshold of 1.3 and 10-second constraint sampling.

To run the commented example, use the existing six-aircraft scenario. Keep vertical_power even so positive and negative altitude differences contribute symmetrically.

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[6],
#     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)