Airport operation robustness decision-making method and device based on stochastic optimization

Abstract

The invention discloses an airport operation robustness decision-making method and device based on stochastic optimization. The method comprises the following steps: constructing a low-altitude intrusion situation feature space through multi-source heterogeneous sensing data; on the basis of the feature space, constructing a runway recovery time Wasserstein uncertainty set capable of dynamically adjusting the boundary so as to describe uncertainty; establishing a distribution robust optimization model for coupling aircraft position distribution and ground support facility scheduling; a column constraint generation algorithm is used for solving, and decision making and closed-loop feedback are executed in a rolling time domain mode. The method does not need to depend on prior probability distribution of interference events, adaptive balance of decision robustness and economical efficiency can be achieved, the problem of space-time mismatching is effectively avoided through resource coupling scheduling, and improvement of the operation recovery efficiency and safety of an airport under high uncertainty is facilitated.

Description

Technical Field

[0001] This application relates to the interdisciplinary fields of artificial intelligence technology, operations research optimization and aviation management, and specifically to a robust decision-making method and apparatus for airport operations based on stochastic optimization. Background Technology

[0002] In the modern airport airspace safety and operations management system, the development of the low-altitude economy and the widespread adoption of unmanned aerial vehicle (UAV) technology have brought new operational risks to airports. Non-cooperative low-altitude flying objects (LAFOs) with high maneuverability, autonomous swarming capabilities, and long endurance may pose an intrusion threat to airport airspace. These sources of interference are fundamentally different from traditional natural weather interference or equipment malfunctions; their behavior patterns are difficult to predict, posing a threat to the normal operation of airports.

[0003] Existing airport operation recovery techniques primarily rely on stochastic optimization theory, such as stochastic programming or robust optimization. On one hand, traditional stochastic programming methods typically assume that the duration of the disruption event or related delays follows a known, fixed probability distribution, such as a normal or Poisson distribution. These methods approximate the optimization problem by statistically analyzing and sampling historical data. However, for non-cooperative low-altitude targets with game-theoretic characteristics, key parameters such as intrusion paths and loiter times exhibit highly non-stationary stochastic process characteristics. The intruder may adjust its strategy in real time based on the airport's countermeasures, causing the probability distribution of relevant random variables to exhibit ambiguity and dynamic time-varying nature. In this situation, the prior distribution fitted based on historical data may become invalid. Continuing to use such erroneous distributions for decision-making will lead to the decision scheme deviating from the optimal state in practical applications, resulting in significant operational costs.

[0004] On the other hand, traditional robust optimization methods construct an uncertainty set encompassing all possible disturbance scenarios and seek the optimal solution under the worst-case scenario within this set to ensure absolute reliability of decisions. However, accurately predicting the maximum potential impact range of intelligent disturbance sources is extremely difficult. To ensure operational safety, the boundaries of the uncertainty set often need to be set quite broadly, such as assuming the runway will be closed for an extended period. This setting can lead to overly conservative decision-making, such as large-scale flight cancellations or prolonged idleness of core resources like runways and aircraft stands, resulting in unnecessary economic losses and a decline in passenger service levels.

[0005] Furthermore, existing technologies suffer from a technical deficiency in decoupling resource allocation when addressing resource scheduling and recovery after runway closure. Specifically, the optimization models for gate allocation and ground support facility scheduling are typically solved independently in stages. This decoupling approach is feasible in deterministic operating environments. However, in high-uncertainty environments caused by low-altitude intrusions, this approach is prone to spatiotemporal mismatches. For example, the command center may allocate available gates to diverted flights based on the gate allocation model, but due to obstructed access or delayed dispatch of ground support vehicles (such as shuttle buses and tow trucks) within the restricted area, the flight may obtain a gate but be unable to receive subsequent ground services, leading to secondary delays and reduced overall operational recovery efficiency.

[0006] Currently, the technology field lacks robust airport operation decision-making capabilities that can effectively overcome dependence on prior distributions and achieve deep coupling and collaborative scheduling of gate resources and ground support resources in scenarios where the probability distribution of interference source behavior is unknown, ambiguous, and dynamically changing. Summary of the Invention

[0007] One aspect of this invention is to provide a robust airport operation decision-making method based on stochastic optimization, which addresses the technical problems of suboptimal decision-making due to reliance on prior distribution and spatiotemporal mismatch caused by resource decoupling scheduling when dealing with interference from non-cooperative low-altitude targets with fuzzy probability distributions.

[0008] To achieve the above objectives, the present invention provides the following technical solution: A robust decision-making method for airport operations based on stochastic optimization includes the following steps: Step 1: Integrate radar data, radio spectrum data, and image data to construct a high-dimensional state vector, which is used to characterize the low-altitude intrusion situation. Step two involves constructing a dynamic Wasserstein uncertainty set for runway recovery time. This step includes: constructing a hybrid density network, using the high-dimensional state vector as input, outputting the weights, mean, and standard deviation parameters of a Gaussian mixture model, and fitting a conditional probability density function for runway recovery time based on these parameters, using the conditional probability density function as a central reference distribution; and quantifying cognitive uncertainty in the prediction using a Bayesian neural network, mapping the prediction confidence interval width, which characterizes the cognitive uncertainty, to the radius of the Wasserstein distance, and dynamically adjusting the boundary of the dynamic Wasserstein uncertainty set using the central reference distribution as the center and the radius as a distance threshold. Step 3: Establish a two-stage sub-Brussels bar optimization model that couples the aircraft position with the ground support facilities. The model includes decision variables for the first stage and decision variables for the second stage, and the objective function is to minimize the sum of the maximum expected values ​​of the deterministic cost of the first stage and the cost of the second stage under all possible probability distributions in the dynamic Wasserstein uncertainty set. Step 4: Solve the two-stage sub-Bruker optimization model of the coupling between the aircraft station and the ground support facilities based on the column constraint generation algorithm. By iteratively solving the main problem and sub-problems, the optimal first-stage decision is obtained. Step 5: Execute the optimal first-stage decision using a rolling time-domain approach, and perform closed-loop feedback recalculation based on real-time situational changes.

[0009] Further, in step one, constructing the high-dimensional state vector includes: defining multiple detected low-altitude targets as nodes in dynamic graph structure data; performing convolution operations on the dynamic graph structure data using a graph attention network to extract topological feature vectors; processing the radar trajectory sequence and spectral fingerprint sequence of the targets using a bidirectional long short-term memory network to extract behavioral pattern feature vectors; and concatenating and fusing the topological feature vectors and the behavioral pattern feature vectors to form the high-dimensional state vector.

[0010] Furthermore, step two also includes: constructing a weighted Mahalanobis distance matrix using the topological feature vector extracted from the graph attention network; defining a weighted Wasserstein distance based on the constructed weighted Mahalanobis distance matrix, thereby constructing a Wasserstein ellipsoid whose shape is determined by the Mahalanobis distance matrix to form the dynamic Wasserstein uncertainty set, wherein the dynamic Wasserstein uncertainty set is an anisotropic uncertainty set.

[0011] Furthermore, the decision variables in the first stage include at least one of issuing instructions to inbound flights to circle or divert to the airport, allocating temporary parking positions to landed flights, and pre-deploying ground support facilities; the decision variables in the second stage include at least one of determining the takeoff order of backlogged flights, changing the parking positions of flights, and planning the service path of ground support facilities.

[0012] Furthermore, when establishing the model in step three, the following constraints are introduced: a service chain time window hard constraint is introduced to ensure that the service start time of ground support facilities is not earlier than the time when the flight arrives at the gate and obtains the gate allocation; and a low-altitude threat-driven dynamic geofencing constraint is introduced to map the threat impact area of ​​low-altitude interference sources to the airport ground road network in real time and restrict the ground support facilities from passing through the threatened area.

[0013] Furthermore, in step four, the main problem is solved by relaxing an infinite number of distributions into a finite set of discrete scenario distributions to obtain the current optimal first-stage decision and the lower bound of the worst-case expected cost; the subproblem, given the first-stage decision, searches for the worst-case probability distribution that maximizes the second-stage expected cost within the dynamic Wasserstein uncertainty set, and obtains the upper bound of the worst-case expected cost.

[0014] Furthermore, solving the subproblem involves using duality theory to equivalently transform the infinite-dimensional optimization problem of finding the worst probability distribution into a finite-dimensional mixed-integer second-order cone programming problem for solution.

[0015] Furthermore, in step five, the execution using the rolling time domain method includes: setting a rolling time window, updating the high-dimensional state vector using the latest sensor data at the end of each time window, and re-executing steps two to four; wherein, the triggering conditions for the closed-loop feedback recalculation include: periodic triggering, i.e., the end of the rolling time window; or event-driven triggering, i.e., detecting a preset situational change in the low-altitude intrusion situation.

[0016] On the other hand, the present invention also provides an airport operation robustness decision-making device based on stochastic optimization, characterized in that it includes: The state vector construction unit is configured to integrate radar data, radio spectrum data, and image data to construct a high-dimensional state vector representing the low-altitude intrusion situation. The dynamic uncertainty set construction unit is configured to construct a hybrid density network, taking the high-dimensional state vector as input, outputting the weights, mean, and standard deviation parameters of the Gaussian mixture model, and fitting a conditional probability density function of the runway recovery time based on the parameters, using the conditional probability density function as the central reference distribution; and using a Bayesian neural network to quantify the cognitive uncertainty in the prediction, mapping the prediction confidence interval width representing the cognitive uncertainty to the radius of the Wasserstein distance, and dynamically adjusting the boundary of the dynamic Wasserstein uncertainty set of the runway recovery time with the central reference distribution as the center and the radius as the distance threshold; The coupling modeling unit is configured to establish a two-stage sub-Brussels bar optimization model that couples the aircraft position with the ground support facilities. The model includes first-stage decision variables and second-stage decision variables, and the objective function is to minimize the sum of the maximum expected values ​​of the deterministic cost of the first stage and the cost of the second stage under all possible probability distributions in the dynamic Wasserstein uncertainty set. The model solving unit is configured to solve the two-stage sub-Bruker optimization model of the coupling between the aircraft station and the ground support facilities based on the column constraint generation algorithm, and obtain the optimal first-stage decision by iteratively solving the main problem and sub-problems. The decision execution and feedback unit is configured to execute the optimal first-stage decision in a rolling time-domain manner and perform closed-loop feedback recalculation based on real-time situation changes.

[0017] The technical solutions provided in the embodiments of the present invention have the following beneficial effects: This invention overcomes the limitations of traditional methods in handling disturbances with unknown, fuzzy, and dynamically changing probability distributions by constructing a data-driven dynamic Wasserstein uncertainty set, achieving an adaptive balance between decision robustness and economy. Simultaneously, by establishing a coupled model of aircraft stands and ground support facilities and introducing dynamic geofencing constraints, it effectively avoids spatiotemporal mismatches in resource scheduling and improves the safety and recovery efficiency of airport ground operations. Finally, by combining an accurate C&CG solution algorithm and a rolling time-domain closed-loop feedback mechanism, a complete robust decision-making scheme for airport operations that can adapt to highly uncertain and adversarial environments is formed. Attached Figure Description

[0018] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0019] Figure 1 This is a flowchart of an airport operation robust decision-making method based on stochastic optimization provided in an embodiment of the present invention.

[0020] Figure 2 This is a schematic diagram illustrating the technology for constructing a low-altitude intrusion situation feature space in an embodiment of the present invention.

[0021] Figure 3 This is a schematic diagram illustrating the principle of constructing a dynamic Wasserstein uncertainty set in an embodiment of the present invention.

[0022] Figure 4 This is a schematic diagram of the iterative solution process of the column constraint generation algorithm in an embodiment of the present invention.

[0023] Figure 5 This is a C-UAS multi-source sensor deployment layout diagram in an embodiment of the present invention.

[0024] Figure 6 This is a spatial structure diagram of low-altitude target cluster cooperative formation in an embodiment of the present invention.

[0025] Figure 7 This is a schematic diagram of the airport ground GSE service network and dynamic geofencing of the present invention.

[0026] Figure 8 This is a comparison chart of the upper and lower bound convergence curves of the C&CG algorithm of this invention. Detailed Implementation

[0027] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0028] Example 1 This embodiment provides a robust airport operation decision-making method based on stochastic optimization. This method can be applied to the decision support system of the Airport Operations Control Center (AOCC) to address highly uncertain events such as intrusions by non-cooperative low-altitude targets. Figure 1 As shown, the specific process of this method is as follows: Step S1: Construct a low-altitude intrusion situation feature space based on multi-source heterogeneous sensor data. The task of this step is to transform real-time, multi-dimensional raw sensor data into structured feature vectors that can comprehensively and accurately characterize the intrusion threat situation. For example... Figure 2 As shown, this step first aggregates data streams from various sensors in real time through the data interface of the counter-drone command system (C-UAS) deployed at the airport. Specifically, data sources may include: S-band or X-band surveillance radars, providing parameters such as target range, azimuth, elevation angle, radius and velocity, and radar cross-section (RCS); distributed radio spectrum monitoring equipment, providing the center frequency, bandwidth, frequency hopping pattern, and received signal strength indication (RSSI) of suspicious communication signals; and electro-optical turrets installed in key areas, providing visible light images and infrared thermal imaging video streams.

[0029] like Figure 5As shown, the deployment layout of the C-UAS multi-source sensors in the airport area in this embodiment of the invention includes the following components: An S / X-band surveillance radar R1 is deployed at a high point in the airport. This radar can provide parameters such as target range, azimuth, elevation angle, radial velocity, and radar cross-section (RCS), and its coverage extends to the entire airport operating area and surrounding airspace. Four radio spectrum monitoring stations F1, F2, F3, and F4 are distributed around the airport to monitor the center frequency, bandwidth, frequency hopping pattern, and received signal strength indication (RSSI) of suspicious communication signals. This distributed deployment enables omnidirectional coverage of the airport perimeter. Electro-optical / infrared turrets E1, E2, and E3 are deployed near both ends of the runway and behind the terminal building, respectively, to provide visible light images and infrared thermal imaging video streams, enabling focused monitoring of the runway approach area and the terminal building area. All sensors are connected to the C-UAS command center via data links to achieve real-time collection and fusion processing of multi-source heterogeneous sensor data. This deployment fully leverages the complementary advantages of the detection characteristics of different types of sensors: radar provides long-range, wide-area detection capabilities, spectrum monitoring equipment provides radio signature recognition capabilities, and optoelectronic turrets provide high-precision target imaging and recognition capabilities. The three work together to achieve comprehensive, multi-dimensional situational awareness of non-cooperative low-altitude targets.

[0030] To effectively characterize potential clustered collaborative attack patterns, this step further abstracts multiple detected suspicious targets into dynamic graph structure data. Specifically, at a certain time t, the detected target i is defined as node vi in ​​the graph. The feature vector of each node can be composed of its kinematic state and physical properties, for example, a seven-dimensional vector containing three-dimensional coordinates (x, y, z), three-dimensional velocity vectors (vx, vy, vz), and RCS values. When the Euclidean distance between any two targets i and j is less than a preset collaborative behavior threshold (e.g., 500 meters), or when a radio spectrum monitoring device detects a clear communication data link between them, an undirected edge e_ij is established between these two nodes vi and vj, forming a dynamically changing graph Gt. Subsequently, a graph attention network (GAT) is used to extract features from graph Gt. Through its built-in self-attention mechanism, GAT can assign different attention coefficients alpha_ij to the neighboring nodes of each node, thus focusing on nodes that are more closely related to the current node's behavior when aggregating neighborhood information. Through stacked computation of multiple layers of GAT, a deep topological feature vector, denoted as h_swarm, is finally generated for the entire cluster, which can characterize its internal topology, formation compactness, motion consistency and potential tactical intentions.

[0031] like Figure 6As shown, this invention abstracts multiple detected suspicious low-altitude targets into dynamic graph structure data to effectively characterize the cooperative behavior of the swarm. In a three-dimensional spatial coordinate system, the X and Y axes represent the horizontal direction, the Z axis represents the altitude direction, and the airport runway is located near the origin as a spatial reference. The runway has entrance and exit markings to enhance identification. The figure shows eight target nodes v1 to v8, each represented by a top-view symbol of a quadcopter UAV, including a central fuselage, cross-shaped arms, and four rotor rings, flying in a typical V-shaped formation. Among them, v1 is the lead node located at the apex of the formation, v2, v4, and v6 are located on the left wing, v3, v5, and v7 are located on the right wing, and v8 is located in the middle rear of the formation, serving as a communication relay. The feature vector of each node includes its three-dimensional spatial coordinates (x, y, z), velocity vector, and radar cross section (RCS) value. An example velocity value of v=25m / s is marked in the figure. Communication data links between nodes are represented by wavy solid lines, reflecting the characteristics of wireless signal transmission; dashed lines represent associations within the cooperative distance threshold but without established direct communication links. The cluster as a whole is surrounded by a dashed elliptical envelope, which indicates the spatial distribution range of the cluster. Velocity vector arrows indicate the direction of movement of each target, suggesting that the cluster is moving towards the airport. The cluster's projection on the ground forms a threat impact area, represented by gray shading and concentric circle warning markers. This area will be mapped onto the airport's ground road network in real time for subsequent dynamic geofencing constraint construction, ensuring the driving safety of ground support vehicles. Feature extraction from this dynamic graph structure using a graph attention network yields deep topological feature vectors characterizing the cluster's degree of cooperation, formation compactness, and potential tactical intentions.

[0032] This step takes the target's radar trajectory sequence (e.g., a sequence of coordinate points over the past 60 seconds) and spectral fingerprint sequence (e.g., a sequence of signal frequency changes over time) as input and feeds them into a bidirectional long short-term memory network (Bi-LSTM). Through its unique gating mechanism, Bi-LSTM can simultaneously learn long-term dependencies in the sequence data from both forward and backward time dimensions, effectively extracting features of behavioral patterns such as hovering, rushing, and evasion, and outputting a behavioral pattern feature vector, denoted as h_behavior.

[0033] Finally, the topological feature vector h_swarm extracted from GAT and the behavioral pattern feature vector h_behavior extracted from Bi-LSTM are concatenated. Further, the concatenated vector is fused with environmental context information obtained from other systems, such as current visibility, wind speed, and wind direction data from a meteorological system, and terminal area flow control level information from an air traffic control system. This information collectively constitutes a comprehensive, high-dimensional intrusion situation state vector St, which serves as input for subsequent steps.

[0034] Step S2 involves constructing a dynamic Wasserstein uncertainty set for runway recovery time. The Wasserstein distance is a distance function used in operations research and statistics to measure the difference between two probability distributions. Based on this distance, the Wasserstein uncertainty set is defined as the set of all probability distributions centered on a reference distribution, covering all distributions whose Wasserstein distance to that reference distribution does not exceed a specific radius. The "dynamic" set constructed in this embodiment is characterized by its center distribution and radius parameters adaptively changing with real-time situational data. The purpose of this step is to construct, without relying on any fixed prior probability distribution assumptions, a data-driven set containing the true probability distribution of runway recovery time tau, i.e., the uncertainty set, using the state vector St. Figure 3 As shown, the construction of this uncertainty set consists of two core steps: the generation of the central reference distribution and the determination of the radius of the dynamic uncertainty set.

[0035] First, a Mixture Density Network (MDN) is constructed to generate the central reference distribution P0. The MDN is a special type of neural network whose output layer does not contain deterministic predicted values, but rather a set of parameters from a Gaussian mixture model. Specifically, the state vector St is used as the input to the MDN, and the network outputs parameters {pi_k, mu_k, sigma_k}, where pi_k is the weight of the k-th Gaussian component, mu_k is the mean, and sigma_k is the standard deviation. Using these parameters, the conditional probability density function of the runway recovery time tau under the current situation St can be constructed, in the form: P0(tau | St) = Σ_k pi_k * N(tau | mu_k, sigma_k^2). Here, N(tau | mu_k, sigma_k^2) represents a normal distribution with mean mu_k and variance sigma_k^2. The distribution P0 integrates historical data and the current situation, and is considered the most probabilistic prediction of runway recovery time, serving as the center of the sub-Bruker optimization model.

[0036] Secondly, to quantify the uncertainty in the current prediction and dynamically adjust the robustness of the decision accordingly, this step utilizes a Bayesian Neural Network (BNN) to assess cognitive uncertainty. Unlike standard neural networks, the weights of a BNN are not fixed values, but rather probability distributions (e.g., Gaussian distributions). During prediction, by sampling the weights multiple times and forward propagating (e.g., using Monte Carlo dropout), the BNN can output a distribution of predicted values ​​rather than a single value. This step uses the BNN to predict the runway recovery time tau and calculates the variance or standard deviation of the output distribution. The magnitude of this variance directly reflects the model's confidence in predicting the current input St: the variance is small when the input St is similar to the training data; the variance is large when St is an unfamiliar situation that the model has never seen before. This step transforms this variance into the radius epsilon_t of the Wasserstein distance using a predefined mapping function (e.g., the linear mapping function epsilon_t = a * var + b). In this way, the radius epsilon_t is dynamically adjusted: when the intrusion behavior pattern is highly predictable, the cognitive uncertainty of the BNN is low, and epsilon_t decreases accordingly, making the subsequent optimization model more focused on economy; conversely, when the intrusion behavior pattern is irregular and unpredictable, the cognitive uncertainty is high, and epsilon_t increases accordingly, the model decision will be more conservative, prioritizing the protection of operational security.

[0037] Furthermore, to more precisely characterize the structure of uncertainty, this step introduces an anisotropic Wasserstein uncertainty set based on cluster topology features. The standard Wasserstein uncertainty set is a sphere centered at P0 with radius epsilon_t, assuming that the true distribution P has an equal probability of deviating from P0 in all directions. However, cluster-induced disturbances may have specific structured risks. Therefore, this step utilizes the cluster topology features h_swarm output by GAT in step S1 to construct a weighted Mahalanobis distance matrix Mt. The construction of this matrix Mt aims to reflect the direction of uncertainty revealed by cluster cooperative behavior. For example, when GAT identifies highly cooperative formations, Mt imposes greater weight on the corresponding distribution moment dimension. Subsequently, this matrix Mt is used to define the weighted Wasserstein distance, thereby constructing a Wasserstein ellipsoid centered at P0, with a shape determined by Mt and a size scaled by epsilon_t. This ellipsoid is stretched in directions of higher risk, providing a greater robustness margin, while relatively contracting in directions of lower risk, avoiding unnecessary conservatism. Finally, the dynamic Wasserstein uncertainty set is defined as: Ft(epsilon_t) = { P ∈ P(Xi) : W_Mt(P, P0)<= epsilon_t}. Here, W_Mt(P, P0) represents the weighted Wasserstein distance between distribution P and the central distribution P0, Xi represents the range of values ​​for the runway recovery time tau, and P(Xi) is the set of all probability distributions on the support set Xi.

[0038] Step S3: Establish a sub-Bruker optimization model that couples aircraft stand allocation with ground support facilities (GSE). This step constructs a model capable of coordinating the optimization of aircraft stand allocation and ground support resource scheduling within the mathematical framework of two-stage sub-Bruker optimization.

[0039] The first-stage decision variable x refers to decisions that need to be determined before the uncertainty (runway recovery time tau) is revealed. These decisions are characterized by pre-deployment and strategic nature. Specifically, x may include: deciding which approaching flights f need to enter the holding area and which need to be diverted to other airports a; allocating temporary parking positions g for flights f that have landed but are stranded on the taxiway due to runway closure; and pre-scheduling critical GSE vehicles v (e.g., towing vehicles, shuttle buses) from their current location to designated standby areas z to shorten subsequent response time. The first-stage cost c'x is the direct cost of these decisions, where ' represents the transpose symbol, such as the additional fuel cost incurred by flight circling, the passenger accommodation and flight repositioning costs resulting from diversions, etc.

[0040] The second-stage decision variable, y(tau), represents the restorative and tactical decisions made based on the specific scenario after the runway recovery time, tau, is observed. These decisions depend on the specific value of tau. Specifically, y(tau) may include: determining the pushback order of departing flights f backlogged due to runway closure; repositioning some flights to optimize subsequent turnaround procedures; and planning specific service routes and schedules for GSE vehicles v for each flight f requiring ground services. The second-stage cost, Q(x, tau), is the cost associated with flight delays, decreased passenger service levels, and GSE operations.

[0041] A key innovation of this step lies in achieving deep coupling between aircraft stand and GSE scheduling by introducing strict constraints. Specifically, for any flight f requiring ground support services (e.g., needing a tow truck to the stand), a hard constraint on the service chain time window is established in the model. This constraint takes the form: T_arrival_f <= T_gate_f <= T_service_start_f. Here, T_arrival_f is the time it takes for flight f to taxi to a waiting point near stand g, T_gate_f is the start time when stand g is officially assigned to flight f, and T_service_start_f is the time when the GSE vehicle providing the service begins operating on flight f. This constraint ensures that a flight must first obtain a stand, and the GSE vehicle must be in place before service can begin, effectively avoiding the spatiotemporal mismatch problem of "stands available but no service".

[0042] To further enhance the model's realism and security, this step also introduces a dynamic geofence constraint driven by low-altitude threats. Specifically, a spatiotemporal network model of airport ground operations is constructed. Arcs in the network represent transfer paths of GSE vehicles between different locations (such as aircraft stands, service points, and standby areas), and the weight of each arc is the transfer time delta_t_ij. When the C-UAS system identifies a threat area of ​​a low-altitude intrusion target, this area is mapped onto the airport ground road network in real time, forming a dynamic geofence. In the spatiotemporal network, any arc crossing this geofence has its passage capacity set to zero or its passage time set to a maximum value during the duration of the threat. This constraint forces GSE vehicles to detour around all threatened areas during route planning, thereby ensuring the physical safety of ground personnel and equipment.

[0043] like Figure 7The diagram illustrates the operational mechanism of the airport's ground-based GSE (Ground-Based Security) service network and dynamic geofencing. The diagram presents the airport's ground layout from a top-down perspective, including the main runway marked with a centerline and runway number (09 / 27), taxiways, and apron areas. The apron area contains multiple aircraft stands (G1-G6), each depicted with an overhead view of the aircraft (wing and fuselage shape), and service nodes (S1-S3). The upper right corner shows the GSE standby area (Z1), which houses ground support facilities such as tow trucks (with booms), shuttle trucks (rectangular with windows and wheels), refueling trucks (with fuel tanks), and baggage carts. Each vehicle is depicted using realistic pictographic icons. Solid lines connect the aircraft stands, service nodes, and standby area, forming the GSE service network. The arcs in this network represent the transfer paths of GSE vehicles between different locations. The dashed circular shaded areas represent low-altitude intrusion threat zones identified by the C-UAS (Continuous Airway Intrusion) system, marked with multiple triangular warning signs to enhance visual alertness. This area is mapped in real time onto the airport's ground road network, forming a dynamic geofence marked "No Entry". When a GSE vehicle needs to travel from standby area Z1 to gate G3 to perform a service task, the section of the original planned route that crosses the threat area is prohibited (indicated by dotted lines and X marks in the diagram). The system automatically plans an alternative route for the GSE vehicle (thick solid arrow in the diagram), detouring from standby area Z1 to the lower right, passing through service node S3 and gate G4 to reach the target gate G3. A towing vehicle icon is drawn on the route to visually indicate the vehicle's driving status. This dynamic geofence constraint mechanism achieves deep coupling between airside recovery decision-making and ground-side physical safety assurance. When a section of road in the threat area is demarcated, its traffic capacity is set to zero, forcing GSE vehicles to detour through the danger zone.

[0044] Finally, the objective function of this bibliometric optimization model can be expressed as: min { c'x + sup E_P [Q(x,tau)]}, with the constraint P ∈ Ft(epsilon_t). Here, ' denotes the transpose sign. The objective aims to minimize the sum of the deterministic cost of the first stage and the maximum value of the expected cost of the second stage under all possible probability distributions P within the uncertainty set Ft.

[0045] Step S4 involves solving the model based on the column-and-constraint generation algorithm. Because the objective function of the model established in step S3 includes an operation of taking the supremum of an infinite-dimensional variable (probability distribution P), it becomes a semi-infinite programming problem and cannot be solved directly. This embodiment uses the column-and-constraint generation (C&CG) algorithm to solve it accurately. Figure 4 As shown, the C&CG algorithm decomposes the original problem into a main problem (MP) and subproblems (SP), and approximates the optimal solution through iterative solving.

[0046] At the start of the algorithm, the main problem (MP) is a relaxed version of the original problem, which only considers a few known worst-case scenarios. Solving MP yields the current optimal first-stage decision x* and the auxiliary variable eta, which is a lower bound estimate of the expected cost in the worst-case scenario.

[0047] The objective of SP is to find the worst probability distribution P* that maximizes the expected cost E_P [Q(x*, tau)] of the second stage, given x*, within the uncertainty set Ft(epsilon_t) defined in step S2. Using duality theory, this infinite-dimensional optimization problem of finding the worst distribution can be equivalently transformed into a finite-dimensional, deterministic mixed-integer second-order cone programming (MISOCP) problem, which can be efficiently solved using commercial optimization solvers.

[0048] If Q*(x*) - eta is less than the preset convergence threshold delta (e.g., 0.1% of eta), it means that the current lower bound eta is close enough to the true upper bound Q*(x*), the algorithm converges, and the current first-stage decision x* is the optimal solution to the original problem. Conversely, if Q*(x*) > eta, it means that the current MP relaxation is too large and needs to be tightened. At this time, the constraint corresponding to the worst distribution P* found from SP (i.e., a Benders cut or an effective inequality) is added back to the constraint set of the main problem. The newly added constraint will cut off a part of the non-optimal solution space, making the lower bound eta obtained in the next iteration monotonically increasing.

[0049] The algorithm iteratively solves the MP and SP problems, continuously adding new "columns" (scenarios) and "constraints" (Benders cuts) to the main problem until the convergence condition is met. Optionally, reinforcement learning (RL) can be introduced as a warm-start strategy to accelerate the convergence of the algorithm. An RL agent is pre-trained offline, which takes the state vector St and the first-stage decision x* as input and quickly outputs the predicted worst-case scenario tau'. When solving the SP problem, this scenario tau' predicted by the RL agent is added as an initial solution or initial scenario, which can guide the search direction of the SP problem, reduce the solution time, and thus reduce the overall number of iterations and computational cost of the C&CG algorithm.

[0050] like Figure 8The figure illustrates the convergence process of the Column Constraint Generation (C&CG) algorithm. The horizontal axis represents the iteration number k, and the vertical axis represents the cost value. The solid line represents the curve of the lower bound, and the dashed line represents the curve of the upper bound Q*. In the standard C&CG algorithm (black curve), the lower bound monotonically increases from an initial low value because the feasible region of the main problem gradually expands with each iteration; the upper bound Q* fluctuates and decreases from an initial high value, reflecting the process of identifying the worst-case scenario in the subproblems. As iterations proceed, the upper and lower bounds gradually converge, and the algorithm terminates when the convergence condition is met, i.e., the difference between the upper and lower bounds is less than a preset convergence threshold. The gray shaded area in the figure represents the convergence region; the standard C&CG algorithm converges on the 15th iteration. Simultaneously, the figure shows the C&CG algorithm accelerated by reinforcement learning (RL) hot-start (gray curve), which reduces the computational cost of solving subproblems by predicting the worst-case scenario, achieving convergence in only 9 iterations, significantly improving algorithm efficiency. This convergence characteristic ensures that the algorithm can obtain the global optimum or its approximate solution within a finite number of iterations.

[0051] Step S5, Rolling Time-Domain Execution and Closed-Loop Feedback. This step transforms the optimization decisions solved by the model into executable operation instructions and establishes a closed-loop feedback mechanism to continuously adapt to the dynamic environment.

[0052] A fixed rolling time window, delta_T, is set, for example, 5 minutes. At the beginning of each time window, the system executes a complete calculation process from step S1 to step S4. At time t=0, the model is solved and the decision within the first time window is executed; at time t=delta_T, the system updates the state vector St+delta_T using the latest sensor data, reassesses the uncertainty and updates the fuzzy set Ft+delta_T, then solves the model again and executes the portion of the new decision for the future delta_T time period. This process continues to roll forward.

[0053] When the C-UAS system detects a significant change in the intrusion situation, such as a sudden increase in the number of targets, a change in cluster formation, or a target entering a pre-defined critical airspace, the system immediately triggers a complete recalculation without waiting for the current time window delta_T to end. This closed-loop feedback mechanism ensures the real-time nature and adaptability of decision-making. For example, if the initial decision is conservative and reserves a lot of resources, but monitoring data shows a decrease in the threat level (epsilon_t decreases) in the next time window, the new optimization decision may release some of the reserved server positions or GSE resources and use them to restore normal operation, thereby improving resource utilization and operational efficiency. Conversely, if the threat suddenly escalates, the decision will also quickly adjust to a more conservative strategy to ensure a safety baseline.

[0054] Example 2 This embodiment provides an airport operation robustness decision-making device based on stochastic optimization. The internal structure of the device can be divided into multiple functional units.

[0055] Specifically, the device includes: The state vector construction unit is used to perform step S1 as described in Example 1. It integrates multi-dimensional sensor data such as radar data, radio spectrum data, and photoelectric or infrared image data in real time through the anti-drone command system. It uses graph attention network to extract cluster topology features, uses bidirectional long short-term memory network to extract behavioral pattern features, and combines meteorological and air traffic control information to jointly construct a high-dimensional state vector that comprehensively represents the current intrusion situation.

[0056] The dynamic uncertainty set construction unit is used to perform step S2 as described in Example 1. Based on the high-dimensional state vector, a central reference distribution is generated using a hybrid density network, and a Bayesian neural network is used to quantify and evaluate cognitive uncertainty to determine the dynamic Wasserstein distance radius. Furthermore, an anisotropic Wasserstein uncertainty set based on cluster topology features is constructed.

[0057] The coupled modeling unit is used to perform step S3 as described in Example 1, and to establish a two-stage sub-Bruker optimization model that couples the allocation of aircraft positions with the scheduling of ground support facilities. The model includes first-stage decision variables and second-stage decision variables, and introduces dynamic geofencing constraints driven by low-altitude threats.

[0058] The model solving unit is used to execute step S4 as described in Example 1, and uses a column constraint generation algorithm to iteratively solve the sub-Bruker optimization model. By alternately solving a main problem and a sub-problem, the algorithm converges, thereby obtaining the optimal first-stage decision.

[0059] The decision execution and feedback unit is used to execute step S5 described in embodiment 1, parse the optimal first-stage decision into business instructions and issue them, and periodically or event-drivenly trigger recalculation in a rolling time domain manner or when a sudden change in the environment is detected, so as to realize closed-loop feedback and dynamic adjustment of the decision.

[0060] Those skilled in the art will understand that the above-described device can be deployed on physical servers, cloud computing platforms, or edge computing devices. Its function is implemented by a processor executing computer program instructions stored in memory. The processor can be a central processing unit (CPU), a graphics processing unit (GPU), or an application-specific integrated circuit (ASIC), etc. The memory can include non-volatile storage media such as random access memory (RAM), read-only memory (ROM), or a hard disk.

[0061] Example 3 This embodiment illustrates the implementation process of the method of the present invention in dealing with low-altitude target swarm intrusion events with high uncertainty and adversarial nature through a specific application scenario, and explains the necessity of adopting the method of the present invention.

[0062] Suppose a large international hub airport's C-UAS system detects eight non-cooperative low-altitude targets flying in formation at the end of the main runway. These targets exhibit erratic circling and maneuvering characteristics. This situation renders traditional stochastic programming methods ineffective because it's impossible to establish a reliable prior probability distribution for the runway closure duration random variable based on historical data. Furthermore, if traditional robust optimization methods are used, to ensure absolute safety, a very long runway closure time (e.g., 2 hours) must be assumed, which would lead to large-scale flight diversions and cancellations, generating unnecessary operational costs.

[0063] The decision-making process of this invention is initiated in this scenario. First, the situation feature space construction unit in step S1 receives multi-source data from the C-UAS system. Radar data shows the kinematic parameters of eight targets, and radio spectrum monitoring confirms the existence of communication links between the targets, indicating that they are a cooperative cluster. The Graph Attention Network (GAT) analyzes the topology of the cluster, and the output topology feature vector h_swarm quantifies its high degree of cooperation. The Bidirectional Long Short-Term Memory Network (Bi-LSTM) analyzes its historical trajectory and identifies its behavior pattern as "persistent, highly maneuverable harassment," rather than simple flight path crossing. These features together constitute the intrusion situation state vector St.

[0064] Next, the dynamic fuzzy set construction unit in step S2 characterizes uncertainty based on the state vector St. The hybrid density network (MDN) generates a central reference distribution P0 based on St with respect to the runway recovery time τ, which may have a mean of 40 minutes. Simultaneously, a Bayesian neural network (BNN) assesses the cognitive uncertainty of the current situation. Since the features extracted by GAT and Bi-LSTM indicate strong unpredictability of target behavior, the BNN outputs a high uncertainty metric, which is then mapped to a large Wasserstein distance radius εt. Furthermore, the specific formation patterns identified by GAT are used to construct an anisotropic Mahalanobis distance matrix Mt, making the final Wasserstein fuzzy set Ft(εt) an ellipsoid, exhibiting greater robustness margin on the risk dimension associated with the formation cooperative attack pattern.

[0065] Subsequently, the decomposed bar optimization modeling unit in step S3 establishes a coupled optimization model. At this point, 15 flights are expected to arrive within the next 30 minutes. The first-stage decision variable x of the model includes issuing circling or diversion instructions to these 15 flights, allocating temporary parking positions to flights that have landed but are blocked, and pre-scheduling ground support vehicles (GSEs). A key coupling constraint is reflected here: the threat area of ​​low-altitude targets is mapped as a dynamic geofence on the ground taxiway network, and any GSE route planning must detour around this area. This avoids the spatiotemporal mismatch problem that may occur in traditional decoupling methods, such as "parking positions have been allocated, but support vehicles cannot reach them due to ground blockages."

[0066] Step S4 initiates the Column Constraint Generation (C&CG) algorithm in the model solving unit. The subproblem searches for the worst-case distribution P* within the dynamic fuzzy set Ft(εt). One potentially identifiable worst-case distribution is a bimodal distribution, meaning the runway has a 40% chance of recovering within 20 minutes (if the target is a feint), but a 60% chance of remaining shut down for 70 minutes (if the target is continuous harassment). The main problem iteratively adds constraints corresponding to this worst-case scenario until convergence, yielding a first-stage decision x* that comprehensively balances operating costs and security risks at the current level of uncertainty.

[0067] Finally, the decision execution and feedback unit in step S5 executes the decision, for example, instructing 10 flights to circle and 5 to divert. The rolling time window ΔT is set to 5 minutes. After 5 minutes, C-UAS detects the cluster target moving towards the airspace east of the airport. The system immediately triggers a recalculation, updating the state vector St+ΔT. As the target's intent and affected area become relatively clear, the cognitive uncertainty of the BNN assessment decreases, leading to a contraction of the fuzzy set radius εt+ΔT. The result of the new round of optimization may be more aggressive, such as canceling divert instructions for some flights and instead continuing to circle and wait, thereby dynamically adjusting the conservatism of the decision and improving the overall recovery efficiency and economy of the airport operation.

[0068] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A robust decision-making method for airport operations based on stochastic optimization, characterized in that, Includes the following steps: Step 1: Integrate radar data, radio spectrum data, and image data to construct a high-dimensional state vector, which is used to characterize the low-altitude intrusion situation. Step two involves constructing a dynamic Wasserstein uncertainty set for runway recovery time. This step includes: constructing a hybrid density network, using the high-dimensional state vector as input, outputting the weights, mean, and standard deviation parameters of a Gaussian mixture model, and fitting a conditional probability density function for runway recovery time based on these parameters, using the conditional probability density function as a central reference distribution; and quantifying cognitive uncertainty in the prediction using a Bayesian neural network, mapping the prediction confidence interval width, which characterizes the cognitive uncertainty, to the radius of the Wasserstein distance, and dynamically adjusting the boundary of the dynamic Wasserstein uncertainty set using the central reference distribution as the center and the radius as a distance threshold. Step 3: Establish a two-stage sub-Brussels bar optimization model that couples the aircraft position with the ground support facilities. The model includes decision variables for the first stage and decision variables for the second stage, and the objective function is to minimize the sum of the maximum expected values ​​of the deterministic cost of the first stage and the cost of the second stage under all possible probability distributions in the dynamic Wasserstein uncertainty set. Step 4: Solve the two-stage sub-Bruker optimization model of the coupling between the aircraft station and the ground support facilities based on the column constraint generation algorithm. By iteratively solving the main problem and sub-problems, the optimal first-stage decision is obtained. Step 5: Execute the optimal first-stage decision using a rolling time-domain approach, and perform closed-loop feedback recalculation based on real-time situational changes.

2. The method according to claim 1, characterized in that, In step one, constructing the high-dimensional state vector includes: Multiple detected low-altitude targets are defined as nodes in the dynamic graph structure data; A graph attention network is used to perform convolution operations on the dynamic graph structure data to extract topological feature vectors; A bidirectional long short-term memory network is used to process the radar trajectory sequence and spectral fingerprint sequence of the target to extract behavioral pattern feature vectors; The topological feature vector and the behavioral pattern feature vector are concatenated and fused to form the high-dimensional state vector.

3. The method according to claim 2, characterized in that, Step two also includes: Using the topological feature vectors extracted from the graph attention network, a weighted Mahalanobis distance matrix is ​​constructed; based on the constructed weighted Mahalanobis distance matrix, a weighted Wasserstein distance is defined, thereby constructing a Wasserstein ellipsoid whose shape is determined by the Mahalanobis distance matrix, to form the dynamic Wasserstein uncertainty set, wherein the dynamic Wasserstein uncertainty set is an anisotropic uncertainty set.

4. The method according to claim 1, characterized in that, The decision variables in the first stage include at least one of issuing instructions to inbound flights to circle or divert flights, allocating temporary parking positions to already landed flights, and pre-deploying ground support facilities; the decision variables in the second stage include at least one of determining the pushback order of backlogged flights, changing parking positions for flights, and planning the service path of ground support facilities.

5. The method according to claim 1, characterized in that, When establishing the model in step three, the following constraints are introduced: A hard constraint on service chain time windows is introduced to ensure that the service start time of ground support facilities is no earlier than the time when the flight arrives at the gate and obtains the gate allocation; In addition, dynamic geofencing constraints driven by low-altitude threats are introduced to map the threat-affected areas of low-altitude interference sources onto the airport ground road network in real time and restrict the access of ground support facilities to threatened areas.

6. The method according to claim 1, characterized in that, In step four, the main problem is solved by relaxing an infinite number of distributions into a finite set of discrete scenario distributions to obtain the current optimal first-stage decision and the lower bound of the worst-case expected cost; the subproblem, given the first-stage decision, searches for the worst-case probability distribution that maximizes the second-stage expected cost within the dynamic Wasserstein uncertainty set and obtains the upper bound of the worst-case expected cost.

7. The method according to claim 6, characterized in that, The solution to the subproblem includes: using duality theory, the infinite-dimensional optimization problem of finding the worst probability distribution is equivalently transformed into a finite-dimensional mixed integer second-order cone programming problem for solution.

8. The method according to claim 1, characterized in that, In step five, the execution using the rolling time domain method includes: setting a rolling time window, updating the high-dimensional state vector using the latest sensor data at the end of each time window, and re-executing steps two to four.

9. The method according to claim 8, characterized in that, The triggering conditions for the closed-loop feedback recalculation include: periodic triggering, i.e., the end of the rolling time window; or event-driven triggering, i.e., a preset situational change is detected in the low-altitude intrusion situation.

10. A robust decision-making device for airport operations based on stochastic optimization, characterized in that, include: The state vector construction unit is configured to integrate radar data, radio spectrum data, and image data to construct a high-dimensional state vector representing the low-altitude intrusion situation. The dynamic uncertainty set construction unit is configured to construct a hybrid density network, taking the high-dimensional state vector as input, outputting the weights, mean, and standard deviation parameters of the Gaussian mixture model, and fitting a conditional probability density function of the runway recovery time based on the parameters, using the conditional probability density function as the central reference distribution; and using a Bayesian neural network to quantify the cognitive uncertainty in the prediction, mapping the prediction confidence interval width representing the cognitive uncertainty to the radius of the Wasserstein distance, and dynamically adjusting the boundary of the dynamic Wasserstein uncertainty set of the runway recovery time with the central reference distribution as the center and the radius as the distance threshold; The coupling modeling unit is configured to establish a two-stage sub-Brussels bar optimization model that couples the aircraft position with the ground support facilities. The model includes first-stage decision variables and second-stage decision variables, and the objective function is to minimize the sum of the maximum expected values ​​of the deterministic cost of the first stage and the cost of the second stage under all possible probability distributions in the dynamic Wasserstein uncertainty set. The model solving unit is configured to solve the two-stage sub-Bruker optimization model of the coupling between the aircraft station and the ground support facilities based on the column constraint generation algorithm, and obtain the optimal first-stage decision by iteratively solving the main problem and sub-problems. The decision execution and feedback unit is configured to execute the optimal first-stage decision in a rolling time-domain manner and perform closed-loop feedback recalculation based on real-time situation changes.

Images