Coordinated Optimization of Inter-Hub eVTOL Feeder Services with Heterogeneous Passenger Behavior

1. Introduction

Electric vertical takeoff and landing aircraft (eVTOL) are increasingly regarded as a promising urban transport option for time-sensitive passengers because they can take off and land vertically, provide point-to-point service, and bypass surface congestion [1,2]. Their deployment also offers a potential route toward more integrated urban transportation networks. Transfers between major transport nodes, such as international airports and high-speed railway hubs, often involve high time costs and strict schedule constraints, making inter-hub feeder service a plausible early application for eVTOL. At the same time, eVTOL deployment must be coordinated with ground transport systems [3], and recent work has begun to integrate UAM into network-equilibrium settings that explicitly account for such multimodal interaction [4]. In a typical air-rail transfer scenario, passengers arriving at a high-speed railway station in the same city may need to reach an airport quickly when direct ground transfer is unavailable or when the transfer window is particularly tight. Related studies on air-to-rail intermodal travel also show that passengers make substantive choices within the transfer chain, including transfer-city selection [5]. In such settings, eVTOL can provide rapid point-to-point feeder service: passengers can transfer from the railway station to a station-integrated vertiport and quickly reach a vertiport retrofitted on the rooftop of an existing parking structure at the destination airport, thereby forming an efficient inter-hub transfer corridor.

Early studies on eVTOL operation and scheduling mainly examined supply-side efficiency under aircraft and infrastructure constraints. Kleinbekman et al. and Pradeep et al. modeled eVTOL arrival sequencing and heuristic scheduling for urban air mobility [6,7]. Other studies incorporated joint routing and charging optimization, time-varying demand, and variable intervals into eVTOL scheduling frameworks [8,9,10]. Farazi and Zou internalized community noise impacts, while Jin et al. considered integrated optimization of strategic planning and service operations for urban air mobility systems [11,12]. These studies primarily focus on supply-side operation or facility planning. By contrast, recent studies in flexible and modular transit systems have started to couple service scheduling with passenger travel-choice preferences or passenger assignment in bilevel frameworks [13,14]. The coupling between passenger choice behavior and operational scheduling therefore remains insufficiently explored.

In practical market operation, operators must consider passenger acceptance of eVTOL. Garrow et al. argued that eVTOL research should jointly consider demand, operations and competition with ground transport [15]. Mode-choice and stated-preference (SP) studies by Boddupalli et al. and Hwang and Hong also show that passenger adoption of eVTOL depends on service attributes, application scenarios and individual preferences [16,17]. More broadly, intermodal and emerging mobility choices exhibit substantial heterogeneity across passenger groups, including air-rail intermodal travelers, MaaS users, and multimodal public transport users [18,19,20]. Jiang et al. demonstrated, using Beijing as an example, that multidimensional demand assessment affects vertiport location [21]. Guo et al. further incorporated uncertain traveler preferences into eVTOL network planning [22]. This motivates a bilevel modeling framework that links operator-side scheduling with passenger-side response. In lower-level demand response, passengers facing a strict remaining connection time constraint and relatively high fares may exhibit highly elastic travel decisions, and travel intentions can shift markedly under flight delays or other disruptions [23]. Their choices depend not only on travel time and monetary cost, but also on direct competition from non-eVTOL feeder alternatives. Evidence from intermodal choice and feeder-mode switching further suggests that such decisions are shaped by a broad set of service and traveler attributes [24,25]. Because of preference differences and perception errors, passenger choice behavior is inherently stochastic. Stochastic user equilibrium (SUE), grounded in utility maximization, introduces perceived error terms and can describe passenger choice and flow distribution more realistically in multimodal transport environments [26]. Incorporating stochastic equilibrium response into eVTOL operation optimization can reduce the mismatch between supply-side schedules and realized passenger allocation, improving both interpretability and applicability.

Although stochastic user equilibrium improves behavioral realism, it substantially increases computational complexity. In multimodal feeder networks, Daganzo and Sheffi identified the structural limitation of the multinomial logit (MNL) model, namely the independence of irrelevant alternatives (IIA) [27]. When alternatives overlap strongly, a simple MNL structure cannot accurately represent correlations among options [28]. In inter-hub feeder settings, correlations often arise both among different eVTOL flights and between eVTOL and ground transport alternatives. If non-additive discrete-choice models such as nested logit are used to represent passenger behavior, standard Karush-Kuhn-Tucker (KKT) conditions cannot easily convert the bilevel model into a tractable single-level problem [29]. Meanwhile, the nested relation between upper-level scheduling and lower-level stochastic equilibrium rapidly increases problem size, limiting the efficiency of traditional exact algorithms and heuristics in large-scale instances.

Data-driven deep neural surrogate models offer a promising way to address this mixed-integer nonlinear bilevel bottleneck [30]. Such methods use deep neural networks to learn implicit mappings in complex optimization models, replacing parts of the analytical structure that are difficult to express explicitly. Dumouchelle et al. proposed the Neural Bilevel Optimization (Neur2BiLO) framework, which approximates the lower-level value function with a neural network and embeds it into upper-level optimization. Fischetti and Jo showed that deep neural networks with rectified linear unit (ReLU) activations can be exactly transformed into mixed-integer linear constraints, providing a methodological basis for embedding surrogate models into optimization problems [31,32]. Machine learning-enhanced Benders decomposition also offers useful ideas for problems that couple large-scale discrete decisions with complex continuous responses [33,34,35]. These algorithmic advances provide a basis for solving large-scale networks while preserving nonlinear lower-level behavioral responses.

This study makes three contributions. First, it develops a passenger eVTOL feeder choice model for inter-hub transfer. Based on an SP survey, the model introduces delay-risk perception under a remaining connection time constraint and combines a mixed nested logit model with a Gaussian mixture model (GMM) to identify heterogeneous passenger groups and their choice responses. Second, it constructs a bilevel eVTOL operation optimization model that accounts for heterogeneous passenger SUE response. The upper-level model maximizes operator profit by selecting eVTOL operation paths and schedules under a given fleet size, charging-infrastructure configuration, and fare level, while the lower-level model describes the equilibrium choices of heterogeneous passengers between eVTOL and external transport modes using nested logit SUE. Third, it proposes a Neur2BiLO-GBD hybrid solution framework and applies it to a representative inter-hub feeder case. By exploiting the separable structure between the linear impedance term and the nonlinear entropy term in the lower-level SUE equivalent objective, the framework shifts the nonlinear entropy term to a relaxed subproblem and iteratively corrects the surrogate master problem through generalized Benders cuts. The case study further evaluates the social net utility of the feeder system after eVTOL entry under profit-maximizing operation and identifies the operating scenarios in which eVTOL offers greater value as a feeder service.

2. Materials and Methods

2.1. Passenger Heterogeneity Modeling

2.1.1. Experimental Design

To examine passenger preferences for eVTOL in inter-hub feeder settings, this study conducted a stated-preference survey built around representative air-rail intermodal transfer scenarios. After arriving at a high-speed railway station in the same city, respondents were asked to choose among feeder modes to reach the airport for an onward flight. The modeled hub transfer process therefore spans arrival at the railway station, feeder-mode choice, and arrival at the destination airport, while subsequent airport processes common to all alternatives are not differentiated.

The core scenario variable is the remaining connection time (RCT), defined as the time available to complete the hub transfer process defined in this study once the passenger arrives at the high-speed railway station and enters the feeder decision stage. The time attribute (Time) of each mode denotes the normal time required to complete that process using the corresponding alternative. The potential delay attribute (Delay) denotes the additional delay that may arise beyond normal operating conditions.

To reduce response burden, the questionnaire directly provided the normal travel time and potential delay of each alternative, rather than asking respondents to calculate the time spent in each transfer stage. Respondents could therefore compare cost, time, and potential delay under a consistent time definition. All passengers were assumed to arrive at the high-speed railway station at the same time, and different remaining connection time levels were used to construct urgency scenarios and identify trade-offs among cost, time, and potential delay. Three scenarios were considered: a relaxed scenario with RCT = 120 min, indicating sufficient time redundancy and a relatively large choice set; a critical scenario with RCT = 90 min, in which the time required under regular feeder options approaches the transfer-time budget; and an urgent scenario with RCT = 75 min, in which the available connection time is further compressed and passengers face higher connection risk.

Because the subsequent optimization stage requires dynamic eVTOL scheduling, different eVTOL flights, as well as eVTOL and ground transport corridors, often overlap in space and time. A nested logit structure is therefore used to capture correlations among alternatives [36]. To represent substitution among service instances within the same mode and to identify the eVTOL nest parameter, two eVTOL alternatives were included in the SP experiment. They represent candidate eVTOL service instances formed by different attribute combinations under a given scenario. The full choice set consists of eVTOL option 1, eVTOL option 2, taxi/ride-hailing, and airport express rail. In the remainder of the paper, "external transport" refers to the non-eVTOL feeder alternatives that are not optimized by the operator, namely taxi/ride-hailing and airport express rail; "ground transport" is used only when describing their physical surface-transport characteristics.

The experimental scenarios were generated in Ngene using a D-efficient design [37]. The final D-error was 0.000124. The SP attributes and levels are listed in Table 1.

The survey targeted passengers in Shanghai with transfer experience involving major integrated transport hubs, such as railway stations and airports. Data were collected through the Tencent Questionnaire platform and its sampling service (https://wj.qq.com/). Before the formal survey, 176 valid pilot questionnaires were collected and used to refine several questionnaire items. The revised questionnaire was then used as the final version. The formal survey was conducted in April 2026. A total of 500 questionnaires were distributed. After screening and validation, 466 valid responses remained. Each respondent completed 18 SP choice scenarios, producing 8388 observations.

Table 2 summarizes respondents' socioeconomic characteristics and psychometric measures. The sample is dominated by passengers aged 26-35 years and by respondents with undergraduate degrees. In terms of travel characteristics, business and official travel account for the larger share. A five-point Likert scale was used to measure psychometric indicators (1=strongly disagree,2=disagree,3=neutral,4=agree,5=strongly agree). The results show high agreement with convenience preference and innovation acceptance, and moderately high levels of price sensitivity, new-technology safety concern, punctuality preference and uncertainty aversion. The sample therefore values efficient direct travel while also paying attention to cost, punctuality and uncertainty risk.

2.1.2. Passenger Choice Model

To estimate passenger preferences across feeder alternatives, this study uses a mixed nested logit model for SP choice data. The model introduces random parameters into a nested logit structure, allowing it to reflect correlations among alternatives and heterogeneity in passenger sensitivities to cost, travel time, and delay-risk perception [36].

According to the four alternatives defined in Section 2.1.1, the two eVTOL alternatives are assigned to the eVTOL nest, while taxi/ride-hailing and airport express rail are assigned to the ground transport nest. Let $m(i)$ denote the nest to which alternative $i$ belongs, and let $B_m$ denote the set of alternatives in nest $m$. Given individual random parameter vector ${\beta }_n$, the probability that passenger $n$ chooses alternative $i$ is:

$$P_{ni}=P_{n,m(i)}{P}_{ni|m\left(i\right)}$$

(1)

$$P_{ni|m}={\displaystyle \frac{\exp \left(V_{ni}/{\lambda }_m\right)}{{\displaystyle \sum _{j\in B_m}\exp }\left(V_{nj}/{\lambda }_m\right)}}$$

(2)

$$P_{nm}={\displaystyle \frac{{\left[{\displaystyle \sum _{j\in B_m}\exp }\left(V_{nj}/{\lambda }_m\right)\right]}^{{\lambda }_m}}{{\displaystyle \sum _q{\left[{\displaystyle \sum _{j\in B_q}\exp }\left(V_{nj}/{\lambda }_q\right)\right]}^{{\lambda }_q}}}}$$

(3)

where $P_{ni|m}$ is the conditional probability that passenger $n$ chooses alternative $i$ within nest $m$, $P_{nm}$ is the probability that passenger $n$ chooses nest $m$, and ${\lambda }_m$ is the inclusive value parameter of nest $m$.

To model delay-risk perception under the remaining connection time constraint, this study introduces a delay-risk perception term in addition to linear time and cost terms [38]. Based on random utility theory [39], the utility of alternative $i$ for passenger $n$ is:

$$V_{ni}=AS{C}_i+{\beta }_{n,c}{C}_{ni}+{\beta }_{n,t}{T}_{ni}+{\beta }_{n,R}{R}_{ni}$$

(4)

where $V_{ni}$ is the systematic utility of alternative $i$ for passenger $n$; $AS{C}_i$ is the alternative-specific constant; $C_{ni}$, $T_{ni}$, and $R_{ni}$ are the cost, travel time, and perceived delay risk associated with alternative $i$ for passenger $n$, respectively; and ${\beta }_{n,c}$, ${\beta }_{n,t}$, and ${\beta }_{n,R}$ are the marginal utility coefficients for cost, travel time, and delay-risk perception. Because these marginal utilities are behaviorally expected to be negative, they are specified as the negative of lognormally distributed random variables.

An exponential delay-risk-perception function is constructed for inter-hub transfer scenarios:

$$R_{ni}=\exp \left(\kappa \cdot {\displaystyle \frac{d_i}{RC{T}_s}}\right)-1$$

(5)

where $d_i$ is the potential delay of alternative $i$, $RC{T}_s$ is the remaining connection time in scenario $s$, and $\kappa$ is the curvature coefficient.

2.1.3. Passenger Heterogeneity Results

To avoid strong correlation between the curvature coefficient $\kappa$ in the delay-risk-perception function and the corresponding marginal utility coefficient ${\beta }_{n,R}$, which may destabilize estimation and create local optima, this study calibrates $\kappa =\{1,2,3,\dots ,16\}$ through grid search and then treats it as a fixed parameter in the main model. As shown in Figure 1(a), model fit is optimized when $\kappa =12$. Model parameters are estimated by maximum likelihood. The full-sample estimation results and goodness-of-fit measures are reported in Table 3.

The main model provides a satisfactory fit and explains mode-choice behavior in this scenario. The mean coefficients on cost, travel time, and delay-risk perception are all significantly negative at the 1% level, indicating that higher feeder cost, longer travel time, and stronger delay-risk perception all reduce passenger utility. The value of ${\rho }^2=0.1723$ is acceptable for discrete-choice models in transport applications.

The standard deviation parameters are also highly significant, indicating substantial unobserved heterogeneity in time valuation and delay-risk perception among passengers facing inter-hub transfer and remaining connection time constraints. To further identify latent behavioral groups, individual posterior preference parameters are extracted and clustered using a Gaussian mixture model (GMM) [40]. Candidate class numbers are set as $G=1,\dots ,6$, and the Bayesian information criterion (BIC) and silhouette coefficient are used for evaluation [41]. Although BIC declines as $G$ increases, the decline shows a clear elbow at $G=2$. The silhouette coefficient also reaches its maximum at $G=2$, indicating the best statistical separation. The final classification therefore uses two classes.

The specific estimates for $G=2$ are shown in Table 4.

The two classes exhibit significant preference heterogeneity under remaining connection time constraints. Their cost coefficients are close, suggesting that, in inter-hub transfer scenarios, the main difference lies in the strength of their sensitivity to travel time and delay-risk perception.

Class 1 accounts for 18.03% of the sample. Its behavioral parameters show stronger negative sensitivity to travel time and delay-risk perception. Its value of time and willingness to pay for unit delay-risk reduction are both higher than those of Class 2. Its $AS{C}_E$ is statistically insignificant, suggesting that, after controlling for cost, travel time, and delay-risk perception, there is insufficient evidence of an additional intrinsic preference for eVTOL. When RCT is reduced from 120 min to 75 min, the eVTOL choice share for this class remains close to 50%. Given its high time value and WTP for delay-risk reduction, this group places greater weight on reducing time loss and the disutility associated with potential delay. We therefore label Class 1 as "robustness-oriented passengers".

Class 2 accounts for 81.97% of the sample. It has lower value of time and lower willingness to pay for delay-risk reduction, but its $AS{C}_E$ is statistically significant, indicating a latent acceptance of eVTOL as a new transport mode. When RCT is reduced from 120 min to 75 min, its eVTOL choice share increases from 35.34% to 75.57%. This class is therefore more sensitive to tightening time constraints and more likely to switch to eVTOL under urgency. We label Class 2 as "general passengers".

2.2. Bilevel Planning Model Considering Heterogeneous Passenger Behavior

The eVTOL operation schedule is formulated as a bilevel network design problem under a leader-follower structure. The operator, as the leader, selects feasible eVTOL operation paths and schedules subject to the given fleet and infrastructure capacities. Passengers, as followers, choose among available modes and flights according to heterogeneous preferences, resulting in a stochastic user equilibrium passenger-flow assignment [42,43]. The resulting passenger allocation directly determines system resource use and revenue.

2.2.1. Physical Scenario

The planning horizon $\mathcal{T}$ is discretized into equal time steps $\Delta t$, where $\Delta t=1$ min. On the supply side, the system operates a homogeneous eVTOL fleet subject to onboard state of charge (SOC), finite take-off/landing resources, and finite charging-pile capacity [8,10,44]. On the demand side, the model represents time-varying bidirectional feeder demand between a pair of hubs. Upon entering the decision stage at the origin hub, each passenger chooses between operator-provided eVTOL service and external transport under a given remaining connection time.

Physical feasibility is enforced through energy continuity, charging duration, and take-off/landing separation. Each flight task $l$ between two hubs is assumed to have fixed flight time and fixed energy consumption $E_l^{flight}$. The remaining energy $E_p(t)$ of the aircraft executing path $p$ at time $t$ must remain within the admissible battery operating range [8,10]:

$$E_{min\text{}}\le E_p\text{}(t)\le E_{max}$$

(6)

In the numerical experiments, $E_{min}=28$ kWh and $E_{max}=140$ kWh, corresponding to a 20% minimum SOC threshold. If the projected energy after the next flight would fall below $E_{min}$, the aircraft must charge before executing that flight [10]. Otherwise, charging can be scheduled as needed. Each charging operation restores the battery to $E_{max}$, and the charging duration $\Delta T_c$ depends on the current remaining energy [8,44]:

$$\Delta T_c=⌈60\cdot {\displaystyle \frac{E_{max}-E_p(t)}{P_{max}}}⌉$$

(7)

where $P_{max}$ is the constant charging power. Because $\Delta T_c$ is measured in minutes and $P_{max}$ is specified in kW, Equation (7) uses the minute-based conversion of the charging time. In the time-space state network, this corresponds to a charging arc from $(h,t,e)$ to $(h,t+\Delta T_c,E_{max})$.

For two consecutive operations $o$ and o' using the same take-off/landing resource, the minimum separation requirement is [10,45]:

$$|t_o^{op}-t_{o^{\prime }}^{op}|\ge T_s$$

(8)

where $t_o^{op}$ and $t_{o^{\prime }}^{op}$ are the times at which operations $o$ and $o^{\prime }$perform take-off or landing.$T_s$ is the minimum separation constant. The binary parameter $u_{h,t,p}^{run}$ indicates whether candidate operation path occupies a take-off/landing resource during the current time step. Given the discrete time step $\Delta t$, the take-off/landing resource capacity of hub $h$, $C_h^{run}$, is:

$$C_h^{run}={\displaystyle \frac{\Delta t}{T_s}}$$

(9)

2.2.2. Time-Space State Network

The time-space state network is defined as a directed acyclic graph $\mathcal{G}=(\mathcal{N},\mathcal{A})$. The node set $\mathcal{N}$ contains nodes $v=(h,t,e)$, where $h\in \mathcal{H}$ denotes a physical hub, $t\in \mathcal{T}$ denotes a discrete time and $e$ denotes battery state of charge. The directed arc set $\mathcal{A}$ contains operation arcs and service arcs.

Operation arcs ${\mathcal{A}}_{phys}$ represent eVTOL circulation and energy transition, including flight arcs, turnaround arcs, charging arcs and parking arcs. A flight arc ${\mathcal{A}}_{fly}$ connects $(h_i,t,e)$ to $(h_j,t+t_l,e-E_l^{flight})$, representing spatial displacement, time progression and energy consumption during a flight task. A turnaround arc ${\mathcal{A}}_{turn}$ connects $(h,t,e)$ to $(h,t+T_{turn},e)$, representing mandatory ground time for boarding, alighting, equipment cooling and safety checks after a flight. A charging arc ${\mathcal{A}}_{chg}$ connects $(h,t,e)$ to $(h,t+\Delta T_c,E_{max})$, representing battery charging at a hub. A parking arc ${\mathcal{A}}_{stay}$ connects $(h,t,e)$ to $(h,t+\Delta t,e)$, representing ground waiting without energy change [10].

Service arcs ${\mathcal{A}}_{serv}$ describe travel choices by heterogeneous passengers in the time-space network, including flight service arcs, external transport arcs, and waiting arcs. A flight service arc ${\mathcal{A}}_F$ connects $(h_i,t)$ to $(h_j,t+t_l)$ and corresponds to a flight arc ${\mathcal{A}}_{fly}$ selected by the upper-level decision. An external transport arc ${\mathcal{A}}_X$ connects $(h_i,t)$ to $(h_j,t+{\tau }_{travel}^m)$, representing completion of the hub transfer by an external mode. A waiting arc ${\mathcal{A}}_{wait}$ connects $(h,t)$ to $(h,t+\Delta t)$, representing time spent waiting inside a hub.

Figure 2 illustrates the operation process of a single eVTOL and the passenger travel-choice process in the time-space network. An eVTOL with initial full energy $E_{max}$ enters the turnaround arc ${\mathcal{A}}_{turn}$ at hub A at time $t_1$, completes ground turnaround in $T_{turn}$, then follows a flight arc ${\mathcal{A}}_{fly}$, flies for $t_l$, consumes $-E_l^{flight}$, and arrives at hub B at time $t_3$. After landing, it enters a charging arc ${\mathcal{A}}_{chg}$, spends $\Delta T_c$ charging and restores its state of charge (SOC) to $E_{max}$.

Passenger choices occur on service arcs ${\mathcal{A}}_{serv}$. Passengers arriving at hub A at $t_1$ choose a feeder mode under the remaining connection time condition $RC{T}_k(t)$. If they choose eVTOL, they first spend waiting time on ${\mathcal{A}}_{wait}$, board at node $(A,t_2)$, and then travel on an eVTOL flight service arc in ${\mathcal{A}}_F$. If they choose external transport, they directly take external transport arc ${\mathcal{A}}_X$ to hub B and complete the hub transfer process. They then proceed through common downstream processes such as security screening and boarding. The labels $x_{a,k}^{r,t*}$ and $c_{a,k}^t$ on eVTOL flight service arcs denote equilibrium passenger flows and perceived impedance solved by the SUE model in Section 2.2.4. Finally, eVTOL passenger flow and external transport flow exit the network at nodes $(B,t_3)$ and $(B,t_5)$, respectively.

2.2.3. Model Assumptions and Notation

Four assumptions are introduced to make the model tractable. First, the eVTOL operation plan is determined uniformly over the planning horizon and does not consider real-time rolling adjustment. Fleet size, take-off/landing resources and charging facilities remain fixed within the planning horizon. Second, passenger demand is aggregated by class. Passengers in the same class share the same utility parameters, and differences across classes are represented by cost, travel-time and delay-risk-perception coefficients. Third, all feeder demand generated during the planning horizon must be assigned to a service. Demand cancellation, trip abandonment and unserved demand are not considered. Fourth, external transport modes do not participate in the operation optimization. Their cost, travel time, waiting time and potential delay are exogenously given. Taxi/ride-hailing and airport express rail waiting times are represented by average waiting times, and external modes are assumed to have unlimited capacity.

Table 5 summarizes the main sets, parameters and decision variables used in the model. Symbols not listed in the table are defined in the text.

2.2.4. Lower-Level Model

The lower-level model is a capacity-constrained stochastic user equilibrium problem. It describes how heterogeneous passengers choose paths in the time-space network according to travel utility under a given upper-level flight supply, and how these choices lead to an equilibrium flow distribution.

The demand density $D_{r,k}^t$ for passenger class $k$ on route $r$ is assumed to follow a two-peak Gaussian model:

$$D_{r,k}^t={\displaystyle \sum _{b\in \{1,2\}}\frac{Q_{r,k,b}}{\sqrt{2\pi }{\sigma }_b}}\exp \left(-{\displaystyle \frac{{(t-{\mu }_b)}^2}{2{\sigma }_b^2}}\right)$$

(10)

where $Q_{r,k,b}$ is the total demand in peak period $b$, and ${\mu }_b$ and ${\sigma }_b$ denote the mean peak time and time coverage span. During hub transfer, each passenger faces a clear remaining connection time. In the dynamic assignment model, $RC{T}_k(t)$ denotes the remaining time budget of passenger class $k$ at decision time $t$, updated from the SP scenario-level $RC{T}_s$ according to the current decision time. This dynamic RCT enters the delay-risk-perception function and affects passengers' perceived utility at each decision time.

The lower-level model uses the utility structure calibrated in Section 2.1 and converts individual preference parameters into class-specific parameters. For passenger class $k$, the utility of choosing eVTOL flight $l$ at time $t$ is:

$$V_{k,l}^E(t)=AS{C}_E+{\beta }_{k,c}{F}_E+{\beta }_{k,t}(t_{arr}^l-t)+{\beta }_{k,R}{R}_{k,l}(t)$$

(11)

where $AS{C}_E$ is the eVTOL alternative-specific constant, $F_E$ is the fixed fare, $t_{arr}^l-t$ is the total travel time including waiting and flight time, and $R_{k,l}(t)$ is the delay-risk-perception function under the passenger's remaining connection time constraint.

According to Equation (5), the delay-risk-perception term is:

$$R_{k,l}(t)=\exp \left(\kappa \cdot {\displaystyle \frac{d_l}{RC{T}_k(t)}}\right)-1$$

(12)

where $d_l$ is the potential delay of eVTOL flight $l$, $RC{T}_k(t)$ is the remaining connection time of passenger class $k$ at decision time $t$, and $\kappa$ is the curvature coefficient.

For external transport mode $m$, the utility function is similarly defined as:

$$V_{k,m}^X(t)=AS{C}_m+{\beta }_{k,c}{F}_m+{\beta }_{k,t}({\tau }_{travel}^m+{\tau }_{wait}^m(t))+{\beta }_{k,R}{R}_{k,m}(t)$$

(13)

where $AS{C}_m$ is the constant of external transport mode $m$, ${\tau }_{travel}^m$ is its base deterministic travel time, ${\tau }_{wait}^m(t)$ is the waiting time at time $t$, and the dynamic delay-risk-perception term is $R_{k,m}(t)=\exp \left(\kappa \cdot {\displaystyle \frac{d_m}{RC{T}_k(t)}}\right)-1$, where $d_m$ is the potential delay of mode $m$.

The generalized utilities are transformed into generalized arc impedances in the time-space graph. For an eVTOL service arc $a$, whose corresponding flight is $l(a)$, impedance is defined as the negative of utility, $c_{a,k}^{E,t}=-V_{k,l(a)}^E(t)$. For an external transport arc $a$, whose corresponding external mode is $m(a)$, impedance is $c_{a,k}^{X,t}=-V_{k,m(a)}^X(t)$.

Based on these generalized arc impedances, the lower-level model uses a nested logit SUE form to describe path assignment by heterogeneous passengers in the time-space service network. Let ${\Omega }_{r,t}$ denote the feasible service-path set for route $r$ and departure time $t$, partitioned into the eVTOL flight nest ${\Omega }_{E,r,t}$ and the ground transport nest ${\Omega }_{X,r,t}$. ${\eta }_k$ is the nest-level scale parameter for passenger class $k$, while ${\mu }_{E,k}$ and ${\mu }_{X,k}$ are the within-nest path-choice scale parameters. In the numerical experiments,${\eta }_k$ is normalized to 1, and the within-nest scales are derived from the full-sample inclusive-value parameters as ${\mu }_{E,k}=1/{\lambda }_e=1.965$ and ${\mu }_{X,k}=1/{\lambda }_g=1$.

Following Sheffi's equivalent mathematical programming approach, under a given upper-level schedule and its resulting flight capacity, the flow distribution of the nested logit SUE model is the optimal solution of the following capacity-constrained convex optimization problem (P1). The probabilistic equilibrium assignment can be expressed through an entropy-optimization form [43,26]. The derivation is given in Appendix A.

$$(P1){\min }_{x,f,y}Z(x,f,y)={\displaystyle \sum _{a\in \mathcal{A}}{\displaystyle \sum _{k\in \mathcal{K}}{\displaystyle \sum _{r,t}{x}_{a,k}^{r,t}}}}{c}_{a,k}^t+{\Psi }_{ent}(f_k^t,y_k^t)$$

(14)

The first term is the generalized total system impedance, representing the cumulative product of arc flow $x_{a,k}^{r,t}$ and corresponding arc impedance $c_{a,k}^t$ over all passengers in the time-space network. ${\Psi }_{ent}$ contains the generalized entropy function with microscopic choice entropy and nested-structure correction:

$${\Psi }_{ent}(f_k,y_k)={\displaystyle \sum _{r,t}[}{\displaystyle \frac{1}{{\mu }_{E,k}}}\mathcal{H}(f_E^t)+{\displaystyle \frac{1}{{\mu }_{X,k}}}\mathcal{H}(f_X^t)+\left({\displaystyle \frac{1}{{\eta }_k}}-{\displaystyle \frac{1}{{\mu }_{E,k}}}\right)\mathcal{H}(y_E^t)+\left({\displaystyle \frac{1}{{\eta }_k}}-{\displaystyle \frac{1}{{\mu }_{X,k}}}\right)\mathcal{H}(y_X^t)]$$

(15)

where $\mathcal{H}(x)={\displaystyle \sum x_j}\ln x_j$ is the entropy operator. $y_{k,r}^{E,t}$ and $y_{k,r}^{X,t}$ are the eVTOL-nest and ground-transport-nest flows for passenger class $k$ on route $r$ at time $t$. Under this scale setting, the entropy coefficients are non-negative and the lower-level equivalent program remains convex.

Let $f_{\omega ,k}^{r,t}$ denote the assignment flow of passenger class $k$ on path $\omega$. Arc flows and nest aggregate flows are obtained from path flows:

$$x_{a,k}^{r,t}={\displaystyle \sum _{\omega \in {\Omega }_{r,t}:a\in \omega }{f}_{\omega ,k}^{r,t}},\forall a,k,r,t$$

(16)

$$y_{k,r}^{E,t}={\displaystyle \sum _{\omega \in {\Omega }_{E,r,t}}{f}_{\omega ,k}^{r,t}},y_{k,r}^{X,t}={\displaystyle \sum _{\omega \in {\Omega }_{X,r,t}}{f}_{\omega ,k}^{r,t}},\forall k,r,t$$

(17)

Given an upper-level schedule $z$, the available seat capacity of eVTOL flight service arc $a$ is determined by the selected operation paths:

$$S_a(z)=S{\displaystyle \sum _{p\in \mathcal{P}}{\xi }_{a,p}}{z}_p,\forall a\in {\mathcal{A}}_F$$

(18)

The model constraints are:

$${\displaystyle \sum _{w:(v,w)\in \mathcal{A}}{x}_{vw,k}^{r,t}}-{\displaystyle \sum _{w:(w,v)\in \mathcal{A}}{x}_{wv,k}^{r,t}}=q_{v,k}^{r,t},\forall v\in \mathcal{N},k\in \mathcal{K},r,t$$

(19)

$$y_{k,r}^{E,t}+y_{k,r}^{X,t}=D_{r,k}^t,\forall k,r,t$$

(20)

$${\displaystyle \sum _{k\in \mathcal{K}}{\displaystyle \sum _{r,t}{x}_{a,k}^{r,t}}}\le S_a(z),\forall a\in {\mathcal{A}}_F$$

(21)

$$x_{a,k}^{r,t}\ge 0,f_{\omega ,k}^{r,t}\ge 0,y_{k,r}^{E,t}\ge 0,y_{k,r}^{X,t}\ge 0$$

(22)

Equations (16) and (17) aggregate path flows into service-arc flows and nest flows, and Equation (18) gives the available capacity of flight service arcs under a given upper-level schedule. Equation (19) is the flow-conservation constraint, where $q_{v,k}^{r,t}$ equals $D_{r,k}^t$ at the origin node, $-D_{r,k}^t$ at the destination node, and 0 at intermediate nodes. Equation (20) assigns all demand to either the eVTOL nest or the ground transport nest, while Equation (21) ensures that passenger flow on eVTOL service arc $a$ does not exceed its available capacity $S_a(z)$. The dual variable ${\gamma }_a\ge 0$ represents the shadow price of capacity scarcity: ${\gamma }_a=0$ when capacity is slack and ${\gamma }_a>0$ when the capacity constraint is binding. Solving the KKT conditions of this convex program, as shown in Appendix A, proves its equivalence to the capacity-constrained nested logit SUE. Equation (22) imposes non-negativity [43].

2.2.5. Upper-Level Model

The upper-level model selects candidate eVTOL operation paths through the binary variables $z_p$, thereby determining the flight schedule under the specified fleet, infrastructure, and fare configuration. The objective is operator profit maximization. $z_p\in \{0,1\}$, $p\in \mathcal{P}$. $z_p=1$ indicates that candidate operation path $p$ is selected, and $z_p=0$ otherwise.

The infrastructure fixed cost $C_{fixed}$, including fleet depreciation and charging-pile construction, is:

$$C_{fixed}=C_{fix}{N}_{ac}+C_{inf}{\displaystyle \sum _{h\in \mathcal{H}}{C}_h^{chg}}$$

(23)

The variable cost $C_{var}(z)$ associated with the operation plan includes flight operating cost and electricity consumption cost. $c_p^{var}$ contains maintenance cost corresponding to flight operating time and the take-off-and-landing cycle cost of each flight leg:

$$c_p^{var}={\displaystyle \sum _{l\in p}(C_{ops}{t}_l+C_{land})}$$

(24)

where $C_{ops}$ is the operating cost per unit flight time, $C_{land}$ is the cost per take-off or landing, and $t_l$ is the duration of a single flight.

$$C_{var}(z)={\displaystyle \sum _{p\in \mathcal{P}}{z}_p}(c_p^{var}+c_{elec}{E}_p)$$

(25)

where $c_{elec}$ is the unit electricity price and $E_p$ is the total energy consumption of path $p$.

Total revenue is determined by the actual flow assigned to each eVTOL flight service arc, $x_{a,k}^{r,t*}$, fed back by the lower-level model, multiplied by the fare. Operator profit $\Pi (z)$ is:

$${\max }_z\Pi ={\displaystyle \sum _{a\in {\mathcal{A}}_F}{\displaystyle \sum _{k\in \mathcal{K}}{\displaystyle \sum _{r,t}{F}_E}}}{x}_{a,k}^{r,t*}(z)-\left(C_{fixed}+C_{var}(z)\right)$$

(26)

The total number of eVTOLs put into operation cannot exceed the operator's fleet size:

$${\displaystyle \sum _{p\in \mathcal{P}}{z}_p}\le N_{ac}$$

(27)

For any hub $h$ and discrete time $t$, the total resource occupation of all selected paths cannot exceed the physical capacity of the hub. The take-off/landing resource and charging-pile constraints are:

$${\displaystyle \sum _{p\in \mathcal{P}}{u}_{h,t,p}^{run}}{z}_p\le C_h^{run},\forall h\in \mathcal{H},t\in \mathcal{T}$$

(28)

$${\displaystyle \sum _{p\in \mathcal{P}}{u}_{h,t,p}^{chg}}{z}_p\le C_h^{chg},\forall h\in \mathcal{H},t\in \mathcal{T}$$

(29)

where $u_{h,t,p}^{run}\in \{0,1\}$ and $u_{h,t,p}^{chg}\in \{0,1\}$ are time-space occupation indicators for take-off/landing and charging resources, respectively. They indicate whether candidate operation path $p$ occupies the corresponding facility at hub $h$ at time $t$. $C_h^{run}$ and $C_h^{chg}$ are the corresponding capacity limits.

2.2.6. Social Net Utility of the Feeder System

Beyond maximizing operator profit at the upper level, this study defines the social net utility of the feeder system to evaluate the combined effect of eVTOL entry on operator profit and improvement in passenger utility. In the no-eVTOL scenario, passengers can only complete inter-hub transfer by external transport modes. For travel by passenger class $k$ on service arc $a$, the evaluation impedance is defined as:

$${\overline{c}}_{a,k}^{eval}=-\left({\beta }_{k,t}{T}_a+{\beta }_{k,c}{F}_a+{\beta }_{k,R}{R}_a\right),a\in {\mathcal{A}}_{serv},k\in \mathcal{K}$$

(30)

where ${\overline{c}}_{a,k}^{eval}$ is the service-arc impedance used for social net utility evaluation, ${\mathcal{A}}_{serv}$ is the set of service arcs, and $F_a$ is the travel cost associated with arc $a$.

Let $G_k^E$ denote the equilibrium evaluation impedance of passenger class $k$ with eVTOL entry, and let $G_k^0$ denote the corresponding impedance without eVTOL entry. The monetized passenger welfare increment is:

$$\Delta C{S}^{eval}={\displaystyle \sum _{k\in \mathcal{K}}\frac{G_k^0-G_k^E}{-{\beta }_{k,c}}}$$

(31)

The social net utility of the feeder system is then the sum of operator profit and monetized passenger welfare increment:

$$\Delta SW=\Pi (z)+\Delta C{S}^{eval}$$

(32)

2.3. Neur2BiLO-GBD Hybrid Solution Framework

Directly solving the eVTOL scheduling problem with stochastic user equilibrium response is difficult for two reasons. First, the problem includes both 0-1 scheduling decisions and nested logit SUE response, creating substantial combinatorial and nonlinear complexity. Second, the lower-level model contains a highly nonlinear equivalent objective with logarithmic entropy terms, so it is not directly amenable to KKT-based linearization. Although generalized Benders decomposition (GBD) is advantageous for high-dimensional discrete optimization, it cannot directly handle nonlinear bilevel nested planning [33,34].

To address these difficulties, this study proposes a generalized Benders decomposition algorithm based on Neur2BiLO, referred to as Neur2BiLO-GBD. The method embeds a neural approximation of the lower-level optimal value function into a single-level surrogate model [31], decomposes that model into a surrogate master problem and a relaxed subproblem, and finally re-solves the original lower-level SUE under the selected schedule to recover exact passenger response.

2.3.1. Single-Level Reformulation Based on Neur2BiLO

The original bilevel model can be written as:

$${\max }_{z,v}\Pi (z,v)=Rev(v)-\left(C_{fixed}+C_{var}(z)\right)$$

(33)

$$\mathrm{s}.\mathrm{t}.z\in \mathcal{Z}$$

(34)

$$v\in \arg {\min }_{v^{\prime }\in {\Omega }_{LL}(z)}Z(v^{\prime })$$

(35)

where $z$ is the upper-level scheduling decision vector, $v=(x,f,y)$ is the lower-level passenger-flow assignment vector, and $\Pi (z,v)$ is operator profit, determined by eVTOL fare revenue $Rev(v)$ and operating cost $C_{fixed}+C_{var}(z)$. $\mathcal{Z}$ is the feasible region of the upper-level physical and logical scheduling constraints. ${\Omega }_{LL}(z)$ is the feasible region of lower-level passenger-flow variables under a given schedule, including flow conservation, demand satisfaction, capacity limits and non-negativity. $Z(v)$ is the equivalent objective of the lower-level nested logit SUE.

Following the value-function reformulation principle proposed by Dumouchelle et al. under the optimistic assumption, the lower-level optimal value function is defined as:

$$\Phi (z)={\min }_{v^{\prime }\in {\Omega }_{LL}(z)}Z(v^{\prime })$$

(36)

Based on this value function, the bilevel model can be rewritten as:

$${\max }_{z\in \mathcal{Z},v\in {\Omega }_{LL}(z)}Rev(v)-\left(C_{fixed}+C_{var}(z)\right)$$

(37)

$$\mathrm{s}.\mathrm{t}.Z(v)\le \Phi (z)$$

(38)

Here, $\Phi (z)$ is the optimal value of the lower-level problem under schedule $z$. Directly computing $\Phi (z)$ requires repeatedly solving the dynamic SUE problem during the combinatorial search, which is computationally expensive. This study therefore uses a multilayer perceptron (MLP) with ReLU activation as the value-function approximation, denoted $\widehat{\Phi }(z;\theta )$, to estimate the lower-level equilibrium value under a given schedule. The surrogate model uses scheduling variables as input and lower-level SUE objective values as training targets. The training samples, network settings, and prediction accuracy are reported in Section 3.2. Following Fischetti and Jo, the forward pass of a ReLU neural network can be represented exactly, via a Big-M formulation, as a set of standard mixed-integer linear constraints with auxiliary continuous and binary variables [31,32].

Because neural-network prediction may contain errors, the surrogate may underestimate the true lower-level value, $\widehat{\Phi }(z)<\Phi (z)$. Since $Z(v)$ is always greater than or equal to the true lower bound $\Phi (z)$, the constraint $Z(v)\le \widehat{\Phi }(z)$ may make the problem infeasible. A non-negative continuous slack variable $s\ge 0$ is therefore introduced, yielding $Z(v)\le \widehat{\Phi }(z)+s$. The upper-level profit objective subtracts a penalty term $\rho s$, encouraging $s$ to approach 0 and forcing the lower-level objective toward the surrogate optimal boundary $\widehat{\Phi }$.

The final reconstructed single-level mixed-integer programming model is:

$${\max }_{z\in \mathcal{Z},v\in {\Omega }_{LL}(z),h,b,\widehat{\Phi },s}Rev(v)-\left(C_{fixed}+C_{var}(z)\right)-\rho s$$

(39)

$$\mathrm{s}.\mathrm{t}.(z,h,b,\widehat{\Phi })\in {\Omega }_{NN}$$

(40)

$$Z(v)\le \widehat{\Phi }+s$$

(41)

$$s\ge 0$$

(42)

where $h$ and $b$ are auxiliary continuous and binary variables derived from neural-network hidden-layer outputs and activation states; ${\Omega }_{NN}$ denotes the linear polyhedral constraints that represent the forward mapping of the surrogate model; $\rho$ is the neural approximation error penalty coefficient; and $s$ is the slack variable. This reformulation converts the original highly nonlinear nested bilevel model into a single-level mathematical programming model with a value-function constraint.

In the studied scenario, solving this single-level model directly with Gurobi remains challenging. First, the constraint $Z(v)\le \widehat{\Phi }(z)+s$ requires the solver to evaluate the lower-level objective $Z(v)$ within the optimization model. Because this function contains high-dimensional continuous flow variables and highly nonlinear logarithmic entropy terms, an MILP solver such as Gurobi would require piecewise-linear or outer-approximation representations for these nonlinear terms. This introduces many auxiliary variables and constraints and causes rapid model expansion. Second, the combination of discrete scheduling decisions and the binary variables introduced by neural-network transformation greatly expands the search space.

2.3.2. Surrogate Master Problem and Relaxed Subproblem

To address the dimensional and nonlinear difficulties of direct single-level solution, this section decomposes the model in Section 2.3.1 into a surrogate master problem and a relaxed subproblem. The master problem retains scheduling decisions, linear passenger-flow relations, resource constraints, and neural surrogate constraints, while using an auxiliary variable to represent the nonlinear entropy term in the lower-level SUE equivalent objective. The relaxed subproblem computes the separated entropy value and gradient at candidate points and returns the corresponding linear cuts to the master problem.

The objective $Z(v)$ in the value-function reformulation constraint can be decomposed into a linear part and a nonlinear entropy term:

$$Z(v)=L(v)+H(v)$$

(43)

where $L(v)$ represents the linear objective term formed by generalized arc impedance and arc flow, and $H(v)$ represents the nonlinear entropy term generated by random utility in the nested logit SUE equivalent form.

The surrogate master problem retains $z$, $v=(x,f,y)$, neural-network linearization variables $h$ and $b$, the value-function prediction $\widehat{\Phi }$, and slack variable $s$. To avoid explicitly expanding $H(v)$ in the master problem, a continuous variable $\alpha$ is introduced to represent the entropy term, and its value is progressively approximated by a set of generated cuts $\mathcal{J}$. The surrogate master problem (MP) is:

$${\max }_{z,v,\widehat{\Phi },\alpha ,s,h,b}Rev(v)-\left(C_{fixed}+C_{var}(z)\right)-\rho s$$

(44)

$$\mathrm{s}.\mathrm{t}.z\in \mathcal{Z}$$

(45)

$$v\in {\Omega }_{LL}^{lin}(z)$$

(46)

$$(z,h,b,\widehat{\Phi })\in {\Omega }_{NN}$$

(47)

$$L(v)+\alpha \le \widehat{\Phi }+s$$

(48)

$$\alpha \ge {\mathcal{C}}^{(j)}(v),\forall j\in \mathcal{J}$$

(49)

$$s\ge 0$$

(50)

where ${\Omega }_{LL}^{lin}(z)$ is the linear feasible region containing flow conservation, path-arc aggregation, demand satisfaction, capacity limits and non-negativity; $\mathcal{J}$ is the current set of generated entropy cuts; and ${\mathcal{C}}^{(j)}(v)$ is the $j$-th entropy cut. Equation (48) is the value-function reformulation constraint in the master problem, and Equation (49) progressively supplements the separated entropy term through linear cuts.

Additional linear strengthening constraints are added to the surrogate master problem to bound serviceable eVTOL demand, flight capacity, and fare revenue, thereby reducing relaxation errors caused by fractional schedules, capacities, and flows.

The valid inequalities added to the master problem correspond to the lower-level demand satisfaction constraint (20), the flight capacity constraint (21) and the upper-level revenue function (26):

$$0\le {\displaystyle \sum _{a\in {\mathcal{A}}_F}{x}_{a,k}^{r,t}}\le D_{r,k}^t,\forall k,r,t$$

(51)

$${\displaystyle \sum _{k\in \mathcal{K}}{\displaystyle \sum _{r,t}{x}_{a,k}^{r,t}}}\le S_a(z),\forall a\in {\mathcal{A}}_F$$

(52)

$$Rev(v)={\displaystyle \sum _{a\in {\mathcal{A}}_F}{\displaystyle \sum _{k\in \mathcal{K}}{\displaystyle \sum _{r,t}{F}_E}}}{x}_{a,k}^{r,t}$$

(53)

where $S_a(z)=S{\displaystyle \sum _{p\in \mathcal{P}}{\xi }_{a,p}}{z}_p$ denotes the available capacity generated by schedule $z$ on service arc $a$.

Let the $k$-th candidate point that triggers relaxed-subproblem checking be:

$$(z^k,v^k,{\widehat{\Phi }}^k,{\alpha }^k,s^k)$$

(54)

where $z^k$ is the candidate schedule, $v^k$ is the candidate passenger-flow assignment, ${\widehat{\Phi }}^k=\widehat{\Phi }(z^k;\theta )$ is the neural value-function prediction, ${\alpha }^k$ is the entropy auxiliary variable and $s^k$ is the slack variable in the value-function reformulation constraint.

The relaxed subproblem fixes candidate passenger flow $v^k$ and exactly computes the nonlinear entropy term separated from the master problem. The entropy value at the current candidate point is:

$$H^k=H(v^k)$$

(55)

The left-hand side of the complete value-function reformulation constraint is:

$$Z^k=L(v^k)+H^k$$

(56)

where $H^k$ is the nonlinear entropy value corresponding to $v^k$, and $Z^k$ is the full value of $Z(v)$ after adding the entropy term back. If the master-problem auxiliary variable underestimates the true entropy term, that is, if $H^k>{\alpha }^k$, a new entropy cut is generated. To construct this cut, the first-order gradient of the entropy term at the current point is computed:

$$g^k=\nabla H(v^k)$$

(57)

Let $v_i$ denote a component of the passenger-flow vector $v$. Then:

$$g_i^k={\left.{\displaystyle \frac{\partial H(v)}{\partial v_i}}\right|}_{v=v^k},i\in {\mathcal{I}}_v$$

(58)

where ${\mathcal{I}}_v$ is the passenger-flow variable index set and $g_i^k$ is the partial derivative of the entropy term with respect to $v_i$ at candidate point $k$. The subproblem outputs $H^k$, $Z^k$ and $g^k$, which are used to determine whether a new cut is needed and to construct it.

2.3.3. Linear Cuts and Algorithm Flow

Using the entropy value $H^k$ and gradient $g^k$ output by the relaxed subproblem, a first-order supporting inequality of the convex entropy function $H(v)$ at $v^k$ can be constructed:

$$H(v)\ge H^k+{(g^k)}^T(v-v^k)$$

(59)

The entropy cut generated at candidate point $k$ is:

$${\mathcal{C}}^{(k)}(v)=H^k+{(g^k)}^T(v-v^k)$$

(60)

When $H^k>{\alpha }^k$, the following linear cut is added to the master problem:

$$\alpha \ge {\mathcal{C}}^{(k)}(v)$$

(61)

Expanding Equation (61) by passenger-flow components gives:

$$\alpha \ge H^k-{\displaystyle \sum _{i\in {\mathcal{I}}_v}{g}_i^k}{v}_i^k+{\displaystyle \sum _{i\in {\mathcal{I}}_v}{g}_i^k}{v}_i$$

(62)

As the cut set $\mathcal{J}$ grows, the auxiliary variable $\alpha$ progressively approaches the entropy term $H(v)$, so that the value-function reformulation constraint in the surrogate master problem becomes increasingly close to that in the original single-level model.

The master-problem decision $z^{*}$ is based on the neural surrogate approximation of the value function and the iteratively generated linear cuts. After the algorithm terminates, the final $z^{*}$ is substituted back into the original lower-level SUE model to compute the true passenger response $v^{*}$ and true profit ${\Pi }_{true}$, thereby correcting the numerical evaluation bias introduced by surrogate approximation and relaxation. The full procedure is shown in Algorithm 1.

Algorithm 1. Neur2BiLO-GBD solution algorithm.

Step	Procedure
Input	Time-space network topology and resource configuration $\mathcal{G}(\mathcal{N},\mathcal{A});$passenger origin-destination (OD) demand set and nested logit parameters; pretrained neural surrogate $\widehat{\Phi }(z;\theta ).$
Output	Scheduling decision $z^{\star };$passenger equilibrium flow assignment $v^{\star };$true system profit ${\Pi }_{true}.$
1	Decompose $Z(v)$as $Z(v)=L(v)+H(v),$and initialize the entropy cut set $\mathcal{J}$
2	Repeat.
3	Construct and solve the surrogate master problem MP, including scheduling constraints, linear passenger-flow relations, neural-network surrogate constraints, linear strengthening constraints and all generated cuts.
4	Obtain candidate point $(z^k,v^k,{\widehat{\Phi }}^k,{\alpha }^k,s^k).$
5	Fix $v^k,$and compute $H^k=H(v^k),$ $Z^k=L(v^k)+H^k$and $g^k=\nabla H(v^k).$
6	If $H^k>{\alpha }^k,$generate entropy cut ${\mathcal{C}}^{(k)}(v)=H^k+{(g^k)}^T(v-v^k)$and add it to the master problem.
7	If the current candidate point satisfies the generated-cut conditions, compute candidate profit ${\Pi }^k=Rev(v^k)-C_{fixed}-C_{var}(z^k)$and update the incumbent feasible solution.
8	Until the master-problem search terminates, the computational budget is exhausted or no effective new cut is generated.
9	Fix $z^{\star }$and solve the original lower-level SUE model to obtain true passenger-flow assignment $v^{\star }.$
10	Compute true system profit ${\Pi }_{true}=Rev(v^{\star })-C_{fixed}-C_{var}(z^{\star }).$
11	Return $(z^{\star },v^{\star },{\Pi }_{true}).$

3. Results and Discussion

3.1. Numerical Experiment Setting

The proposed model and solution method are evaluated using a representative air-rail intermodal feeder scenario between two major hubs. The experiments first compare Neur2BiLO-GBD, Neur2BiLO, mixed-integer Karush-Kuhn-Tucker reformulation (MKKT) and a genetic algorithm (GA), which is commonly used for bilevel planning models, under different problem scales. This comparison evaluates computational efficiency and solution quality. The baseline eVTOL schedule and passenger-flow allocation are then analyzed to illustrate system operation under complex supply-demand conditions. The effects of fleet size and charging-facility configuration on operator profit and social net utility are examined. Finally, operational outcomes under different fares, remaining connection times, and external transport delay levels are analyzed.

The case study uses the feeder corridor between Shanghai Hongqiao Integrated Transportation Hub and Pudong International Airport. The straight-line distance is approximately 45 km, and the actual ground feeder distance is approximately 60 km. This corridor is representative of long-distance dual-hub transfer in a large city. Passenger flow at major transport hubs is strongly concentrated in time and space, making it suitable for inter-hub feeder service [46]. Because current eVTOL operations are constrained by visual flight rules (VFR) and visual meteorological conditions, safe night operations under low visibility are difficult to ensure. The eVTOL operating window is therefore set to 08:00-18:00 [47].

Figure 3. Inter-hub feeder scenario between Hongqiao Hub and Pudong International Airport. 

Figure 3. Inter-hub feeder scenario between Hongqiao Hub and Pudong International Airport. 

Experimental parameters include operator supply-side parameters, external transport parameters, and passenger behavior parameters. On the supply side, a representative five-seat eVTOL is used to represent a homogeneous fleet. Single-flight duration is set to 14 min, and single-flight energy consumption is set to 30 kWh. Other parameters are given in Table 6. Passenger behavior parameters are taken from the estimation results in Section 2.1. Passengers are divided into robustness-oriented passengers and general passengers, with market shares of 18.03% and 81.97%, respectively. The cost, travel-time, and delay-risk-perception coefficients of each class follow Table 4, and the curvature coefficient in the delay-risk-perception function is set to $\kappa =12$. Total daily passenger demand is allocated to the two classes according to these market shares.

On the demand side, the alternative set is constructed from the travel characteristics of Hongqiao-Pudong transfer passengers. Existing air-rail intermodal studies show that transfer time, feeder convenience, intermodal passenger behavior, and timetable coordination all affect air-rail connection service quality [48,49,50]. Such passengers usually have high values of time, face tight schedule constraints, and may carry luggage, which makes them particularly sensitive to feeder time and travel reliability. This study therefore excludes conventional urban rail transit from the external alternatives because it is slower, involves more stops, and offers lower comfort. Instead, taxi/ride-hailing and airport express rail are selected as the main external substitutes for eVTOL. Taxi/ride-hailing provides door-to-door feeder service, whereas airport express rail provides a rapid rail feeder connection between Hongqiao Hub and Pudong Airport.

Service and cost parameters for each alternative are set using average road conditions from AutoNavi, historical fare data, and public operating information. In the baseline scenario, taxi/ride-hailing is assigned an average travel time of 60 min through the urban expressway network, an average waiting time of 5 min, and a fare of 150 CNY. Airport express rail is assigned a travel time of 40 min, an average waiting time of 7.5 min, and a full-trip fare of 26 CNY. These time parameters correspond to the hub transfer process defined in Section 2.1 and exclude common subsequent airport processes. The baseline configuration sets the eVTOL fleet size to 16 aircraft, the total number of charging piles in the two-hub system to 8, total daily passenger demand to 4000 passenger trips, remaining connection time (RCT) to 100 min, and eVTOL fare to 200 CNY.

3.2. Algorithm Performance

To evaluate solution performance, this section uses MKKT as the baseline algorithm and compares Neur2BiLO, Neur2BiLO-GBD and GA across different computational scales. The eVTOL operating durations are 60 min, 120 min, 240 min, 360 min and 600 min. All methods use the same baseline configuration, and the time limit for each instance is 600 s. All numerical experiments are run on a computer with a 13th Gen Intel Core i7-13700HX central processing unit (CPU) at 2.10 GHz and 32 GB random-access memory (RAM).

For Neur2BiLO and Neur2BiLO-GBD, independent MLP surrogate models are trained for different-scale instances to approximate the lower-level value function. Training samples consist of 30,000 feasible schedules and their corresponding lower-level SUE solutions. Scheduling decision variables are used as model inputs, and lower-level value-function values are used as fitting targets. Samples are split into training, validation and test sets at a ratio of 8:1:1. Each surrogate model has two hidden layers with 96 neurons per layer, uses ReLU activation, and is trained with the Adam optimizer at a learning rate of 0.01. Test results show high predictive accuracy for all scales, with test-set coefficients of determination above 0.98. For the 600-min instance, the test-set $R^2$ is 0.9887, and the mean absolute error (MAE) and root mean square error (RMSE) are 232.9 and 350.0, respectively.

As shown in Table 7, MKKT obtains the optimal solution in the 60-min case. As operating duration increases, however, the complementarity constraints and integer variables introduced by KKT single-level reformulation cause the model size to grow rapidly, and MKKT cannot find feasible solutions within the time limit for larger cases. GA can return feasible schedules within the time limit, but each individual evaluation requires a passenger-flow assignment procedure. As the search dimension grows with operating duration, the effective search depth becomes limited. Its solution quality therefore fluctuates, and it cannot provide an optimality gap.

Neur2BiLO rapidly obtains objective values identical or very close to MKKT in the 60-min and 120-min cases, indicating that value-function-based single-level approximation can reduce bilevel solution difficulty. As duration grows, however, piecewise linearization of the SUE entropy term introduces many additional variables and constraints, and solution efficiency declines. In the 240-min and larger cases, Neur2BiLO shows longer solution times and larger optimality gaps, suggesting that direct solution of the single-level reformulation remains computationally demanding at large scale. In the largest instance, the large gap of direct Neur2BiLO mainly reflects loose bounds caused by the expanded single-level reformulation, whereas Neur2BiLO-GBD maintains a much tighter bound. By separating the SUE entropy term into a subproblem and using Benders cuts to correct the master problem, Neur2BiLO-GBD avoids most of that expansion and remains more stable in larger cases. For the 240-min, 360-min, and 600-min cases, it obtains higher-quality feasible solutions and smaller optimality gaps.

Figure 4 shows the normalized upper and lower bound convergence trajectories of Neur2BiLO-GBD for representative 60-min and 240-min instances. In the 60-min case, the bounds converge quickly and approach the optimum within a short time. In the 240-min case, the bounds gradually tighten as Benders cuts are added, indicating that the decomposition algorithm effectively uses subproblem information to correct the master problem. For the 600-min large-scale instance, the GBD process must iteratively generate cuts to correct the master problem, and the larger operating scale makes the master relaxation harder to approximate. The optimality gap therefore cannot fully close within the time limit. Nevertheless, Neur2BiLO-GBD stably obtains high-quality feasible solutions and provides valid optimality bounds under the given computational budget.

3.3. Scheduling Results

To illustrate baseline scheduling outcomes, this section examines flight execution, charging-pile occupation, passenger service, and representative aircraft trajectories.

Figure 5 and Figure 6 show that the baseline timetable concentrates flight supply in the morning and evening demand peaks while preserving off-peak windows for charging and operational adjustment. No take-off or landing operation violates hub resource-capacity constraints, and charging-pile occupation remains within the system-wide capacity of the eight charging piles. The resulting schedule therefore maintains fleet circulation while satisfying the physical resource constraints.

Passenger arrivals in Figure 7 exhibit a clear morning-evening double-peak pattern. eVTOL ridership rises with those peaks, but fleet size and seat capacity still leave a substantial share of demand to external transport. In total, 1183 passengers choose eVTOL, and Class 2 passengers account for the larger share of eVTOL ridership.

Figure 8 illustrates this coordination at the aircraft level. For AC09 and AC10, repeated flight tasks are sustained by multiple off-peak charging events, with SOC always remaining within the admissible range and without prolonged ineffective idle periods. These trajectories are consistent with a feasible timetable that coordinates demand peaks, charging, and fleet turnaround.

3.4. Resource-Configuration Effects

Under baseline fare, passenger demand, and external transport parameters, Figure 9 reveals a capacity-matching effect between fleet size and charging infrastructure. Operator profit increases as fleet size rises from 12 to 16 aircraft, but declines when the fleet is expanded to 24 aircraft. For a fixed fleet size, adding charging piles is also beneficial only up to a point: with 16, 20, and 24 aircraft, profit rises when charging piles increase from 6 to 8, then falls as further charging investment is added; with 12 aircraft, profit declines as charging-pile capacity increases. Moderate increases in flight and charging capacity therefore improve service supply, whereas excessive resource investment reduces utilization and raises fixed costs.

The social net utility of the feeder system remains positive under all tested configurations, indicating that eVTOL entry yields net positive benefits across the scenarios examined here. Compared with operator profit, social net utility is more stable at larger fleet sizes. The 24-aircraft configuration yields lower profit than the 16-20-aircraft cases, yet still maintains relatively high passenger-utility gains. Comparing the two heatmaps, the 20-aircraft, 8-pile configuration lies in a high-profit region and yields the highest social net utility in this experiment. Under the Hongqiao-Pudong setting, this combination best balances operating return and passenger-utility improvement, whereas the 24-aircraft configuration still improves passenger utility but reduces social net utility because of higher fixed investment.

3.5. Operational Scenario Analysis

This section analyzes the effects of fare and remaining connection time under the baseline fleet and charging-facility configuration. Figure 10 shows that these two factors shape eVTOL feeder-service performance through different mechanisms. As fare increases, operator profit first rises and then falls, with a high-profit region near 250 CNY, whereas social net utility remains high around 150-200 CNY before declining. From a feeder-system perspective, fare changes reallocate the value of eVTOL service between operator revenue and passenger utility: higher fares increase revenue per passenger but also raise generalized travel cost and weaken passenger-utility gains. The fare that maximizes profit therefore does not coincide with the fare that maximizes social net utility. The 100 CNY scenario does not yield the highest social net utility either, because excessively low fares weaken operator revenue and reduce the overall system improvement. When fare is further increased to 300 CNY, operator profit and social net utility both decline because the demand base contracts.

Remaining connection time affects how fully eVTOL's time and reliability advantages can be translated into realized value. When RCT increases from 80 min to 120 min, both operator profit and social net utility decline markedly. Under a tighter remaining connection time, passengers face stronger time pressure during the inter-hub transfer, so eVTOL attracts more demand through its travel-time and reliability advantages. As remaining connection time relaxes, taxi/ride-hailing and airport express rail can satisfy more feeder demand, and the relative advantage of eVTOL weakens. In the current case, moderate fares and tighter transfer-time constraints are more favorable to higher social net utility of the feeder system.

To further examine how external transport reliability affects eVTOL feeder-service performance, this section extends the baseline setting by varying remaining connection time, fare, and the potential delay of external transport. As defined in Section 2.1, longer potential delay and shorter RCT strengthen delay-risk perception. The potential delay of eVTOL is fixed at 5 min. In the low-delay scenario, the potential delay of airport express rail and taxi is set to 5 min and 10 min, respectively; in the high-delay scenario, it is increased to 10 min and 20 min.

Higher external delay strengthens the value of reliable aerial feeder service, especially when RCT is short. Under low delay, higher-profit outcomes are concentrated near 200 CNY. Under higher external delay, the 200-300 CNY fare range maintains relatively high profit at RCT = 80 min and RCT = 100 min. Social net utility follows a different pattern. Under low delay, the high-value region is concentrated mainly in the 100-200 CNY range. Under high delay, even the 300 CNY scheme remains in a region of relatively high social net utility at RCT = 80 min. At RCT = 100 min, the fare range associated with relatively high social net utility is concentrated near 200 CNY, whereas at RCT = 120 min, the 100 CNY scheme yields the highest social net utility. As remaining connection time relaxes, the high-value fare range shifts back toward lower fares.

This pattern arises because passengers place greater value on reliable feeder service when remaining connection time is short, allowing the time savings and reliability gains offered by eVTOL to sustain both higher fares and high social net utility. As remaining connection time relaxes, the marginal value of reliability improvement declines. Profit peaks therefore remain concentrated at medium-to-high fares, but the fare range associated with high social net utility shifts downward.

Figure 11. eVTOL feeder-service performance under external-delay scenarios.

Figure 11. eVTOL feeder-service performance under external-delay scenarios.

4. Conclusions

This study develops a bilevel planning model for eVTOL feeder operations in transfer scenarios between integrated transport hubs and proposes a Neur2BiLO-GBD hybrid solution framework to improve solution efficiency. The upper level selects eVTOL operation paths and schedules under specified fleet, charging-facility, and fare configurations, while the lower level captures heterogeneous passenger choices among feeder alternatives.

First, passenger choice of eVTOL feeder service shows clear heterogeneity. The parameter estimation results indicate that heterogeneity is mainly reflected in sensitivity to travel time and delay-risk perception, while differences in price sensitivity are relatively small. Numerical experiments further show that fare changes mainly affect total eVTOL demand, whereas transfer-time conditions and external transport reliability more strongly amplify differences in passenger response. Considering passenger heterogeneity in inter-hub feeder problems therefore provides a more realistic representation of lower-level mode choice and helps identify the transfer scenarios in which eVTOL offers greater value as a feeder service.

Second, the proposed Neur2BiLO-GBD framework improves the solution efficiency of the bilevel optimization model. Numerical experiments show that, in small-scale cases, the method obtains results identical or very close to those of MKKT. As instance size increases, it delivers higher-quality feasible solutions within limited computation time and provides corresponding optimality bounds compared with direct single-level reformulation and GA. Combining neural value-function approximation with Benders decomposition therefore provides an effective solution approach for bilevel optimization problems that contain passenger equilibrium response and integer scheduling decisions.

Third, in the Hongqiao-Pudong inter-hub feeder case, eVTOL operational performance is strongly scenario-dependent. Increasing fleet size and charging facilities improves service capability, but marginal returns decline once supply exceeds actual demand. Fare level, remaining connection time, and external transport conditions jointly affect eVTOL market attractiveness. When passengers face tighter transfer-time constraints or lower reliability in external transport, the timeliness and reliability advantages of eVTOL are more likely to translate into higher operator profit and greater social net utility for the feeder system. In low-pressure transfer scenarios, external transport can satisfy most feeder demand, and eVTOL mainly plays a supplementary role. eVTOL therefore offers greater value as a feeder service in inter-hub scenarios with high time sensitivity and high reliability demand.

Fourth, operator profit and social net utility can differ substantially under different pricing schemes. The external-delay scenario analysis shows that, under high delay, regions of higher operator profit are concentrated in medium-to-high fare ranges, while the fare range with high social net utility shifts toward lower fares as remaining connection time becomes more relaxed. Operator profit mainly reflects the financial return from reliable feeder service, whereas the social net utility of the feeder system is jointly determined by operator profit and changes in passenger utility. When transfer time constraints are strong, medium-to-high fares can achieve both high operator profit and high social net utility. As remaining connection time relaxes, however, the fare range with high social net utility gradually falls below the fare range associated with high operator profit. eVTOL operating strategies should therefore consider both operator profit and social net utility, with fare and resource configurations selected according to the specific operating objective.

This study has several limitations. First, the numerical experiments are based on a single inter-hub feeder scenario and do not consider multi-origin-destination (OD) demand, node coordination or route-network reconstruction in a multi-hub urban agglomeration network. Second, weather disturbances, airspace capacity, route conflicts and real-time passenger arrival fluctuations are represented in a simplified manner. Third, passenger behavior parameters are mainly obtained from SP survey data and should be updated and calibrated using revealed-preference data once real operations become available. Fourth, from an algorithmic perspective, the current decomposition framework still has weak mixed-integer relaxation bounds in large-scale instances. Because 0-1 decisions and passenger-flow assignment are tightly coupled, the root-node continuous relaxation cannot sufficiently approximate the integer feasible region, making the optimality gap difficult to close within limited time. Future work can explore strengthened formulations based on path columns, aircraft-rotation columns or time-space arc flows, and combine them with column generation, branch-and-price or logic-based Benders decomposition to improve relaxation tightness and solution stability in large-scale eVTOL feeder networks.