Integrating Flight Scheduling, Fleet Assignment, and Aircraft Routing Problems with Codesharing Agreements under Stochastic Environment

Abstract

Airlines face the imperative of resource management to curtail costs, necessitating the solution of several optimization problems such as flight planning, fleet assignment, aircraft routing, and crew scheduling. These problems present some challenges. The first pertains to the common practice of addressing these problems independently, potentially leading to locally optimal outcomes due to their interconnected nature. The second challenge lies in the inherent uncertainty associated with parameters like demand and non-cruise time. On the other hand, airlines can employ a strategy known as codesharing, wherein they operate shared flights, in order to minimize these challenges. In this study, we introduce a novel mathematical model designed to optimize flight planning, fleet assignment, and aircraft routing decisions concurrently, while accommodating for codesharing. This model is formulated as a three-stage non-linear mixed-integer problem, with stochastic parameters representing the demand and non-cruise time. For smaller-scale problems, optimization software can effectively solve the model. However, as the number of flights increases, conventional software becomes inadequate. Moreover, considering a wide array of scenarios for stochastic parameters leads to more robust results; however, it is not enabled because of the limitations of optimization software. In this work, we introduce two new simulation-based metaheuristic algorithms for solving large-dimensional problems, collectively called “simheuristic.” These algorithms integrate the Monte Carlo simulation technique into Simulated Annealing and Cuckoo Search. We have applied these simheuristic algorithms to various problem samples of different flight sizes and scenarios. The results demonstrate the efficacy of our proposed modeling and solution approaches in efficiently addressing flight scheduling, fleet assignment, and aircraft routing problems within acceptable timeframes.

1. Introduction

The airline industry is highly competitive, especially with the emergence of low-fare airlines. It is characterized by variable demand, high operating costs, heavy traffic, and strict regulations, all of which contribute to the tight competition in the market. Furthermore, an airline’s sole product is its aircraft seats, which are not sold until the flight takes off. Therefore, managing supply and demand fluctuations effectively is a significant challenge for airline companies. Accurate models and solution methodologies with narrow margins of error are necessary to address this issue [1].

Airline operation planning is a complex process of making decisions with many important decisions that directly affect a company’s profitability. The planning process is divided into sub-problems to manage these decisions [2]. The flight scheduling (FS) problem defines the frequency of flights to each destination and how the flights should be planned to meet this frequency. The fleet assignment (FA) problem assigns fleets to appropriate flights, while the aircraft routing (AR) problem creates flight route plans for each aircraft, considering that maintenance constraints. The crew scheduling (CS) problem involves determining the crew members for each flight [3].

Airline companies aim to expand their airline networks to provide passengers with a wider range of destinations, regional connections, and seamless travel experiences with minimal schedule disruptions. However, achieving these goals often involves large investments in purchasing or leasing new aircraft and hiring additional crew members. Alternatively, airline companies can establish codeshare agreements with each other to reduce costs and improve their network coverage.

A codeshare agreement (CA) is a business arrangement between two airlines. One of them can market a flight that is operated by another airline on or off its network, without using their own aircraft. This airline is called a marketing carrier. The other airline operating the flight is called the operating carrier [4,5]. These agreements are made as a result of alliances between two companies, and they allow airlines to expand their flight targets without incurring the huge costs associated with acquiring new aircraft. CAs provide more flight opportunities within an airline’s network, which benefits both customer satisfaction and reduces delays. Additionally, such agreements lead to increased mutual revenues and market shares for the companies [6].

This study aims to simultaneously solve the FS, FA, and AR problems considering CA. It also creates airline schedules that affect profitability while assigning the most suitable fleet types to flights. This improves the company’s service quality and profitability.

The paper is organized as follows: Section 2 includes a literature review of the FA, FS, and AR problem. Section 3 defines the problem and provides a new mathematical model. Section 4 introduces the proposed solution methods and includes a small sample. Section 5 presents computational experiments. Finally, Section 6 concludes the article and discusses future studies.

2. Literature Review

Airlines aim to maximize their profits by providing the optimum combination of fleet type and passengers’ demand. Assigning an aircraft that is smaller than the required passenger capacity to a flight leads to a direct loss of customers due to capacity insufficiency. Conversely, assigning an aircraft that is larger than the required passenger capacity to a flight can result in unsold seats. As a result of this, it causes an inability to sell the seats and reduce operating profits. Therefore, the FA problem is a crucial decision for an airline’s planning strategy [7]. With this motivation, this study examines the FA problem, and Table 1 summarizes the features that contribute to the development of the problem and the studies that investigate these features.

There have been numerous studies on the FA problem in the literature since 1971. While the models of Abara [8] and Hane et al. [9] have contributed to the development of solutions to the problem, their applicability is limited due to certain assumptions such as fixed departure times for the same flight every day and deterministic passenger demand, which do not reflect real-life variability. Clarke et al. [10] developed a solution by incorporating the maintenance planning (MP) of aircraft at stations and crew issues into the FA model.

Berge and Hopperstad [11] suggested the re-flying concept to address demand fluctuations and used the demand-based dispatch approach to make capacity flexible allowing for fleet type assignments to occur closer to take-off. Talluri [12] extended their work in terms of opportunities and computation times. Levin [13] provided flexibility to the FA problem by allowing for variable take-off times and he was the first researcher to consider varying departure times with binary constraints by using integer programming. However, this model did not exclude the aircraft capacity and multiple fleet types. Rexing et al. [14] increased the flight connectivity possibilities by using a time window (TW) for departure times and optimizing time windows according to specific time zones. Belanger et al. [15] minimized the number of aircraft in the fleet by determining the fleet size (FSD). Thus, it aims for the utilization of the fleet (FU) optimally. ZEGHAL ET AL. [16] have addressed a flexible aircraft fleeting and routing problem, which is motivated by the Tunisian national carrier TunisAir. Flexible fleet (FF) is the ability to expand the fleet by renting an aircraft. Thus, more flights can be served.

Barnhart et al. [17] considered the passenger demand based on an itinerary (IFAP) and passenger spill and rescue costs and analyzed the profitability of the model by taking into account the effect of the flight network. In real-life scenarios, unexpected situations and disruptions can occur, so it is necessary to create robust models. Rosenberger et al. [18] developed a robust model using the hub-and-spoke network structure and produced short cycles sensitive to flight cancellations. The robustness of such an assignment and rotation has been demonstrated by using simulations of airline operations. Cadarso and Marín [19] proposed a robust model for the integrated FS and FA problem that takes into account the number of misconnected passengers. The model uses the exponential distribution of passenger connection times (PCTs) to calculate the probability of passengers missing their flight.

Smith and Johnson [20] proposed the concept of station purity (SP) by limiting the number of fleet type or crew-compatible families. Because each airport cannot serve each fleet type, Sherali et al. [21] developed a demand-driven re-flight model that dynamically reassigns the fleet type. Jacobs et al. [22], Dumas et al. [23], and Pilla et al. [24] relaxed the deterministic assumptions set by the basic FA problem. Sherali and Zhu [25] developed a two-stage stochastic programming model for the FA problem. Pilla et al. [26] explicitly evaluated the passenger demand by including direct demand scenarios in Sherali and Zhu’s [25] model. Naumann et al. [27] explained the uncertainty in demand as well as the uncertainty in fuel prices. Cadarso and Celis [28] developed a stochastic model that considers demand uncertainty (SD) and discrete passenger choice models for the FS and FA problem. The model aims to maximize the airline’s profits while minimizing the number of misconnected passengers. The passenger spill and rescue effects are also modeled and it focuses on the probability of passengers choosing a route. Atasoy et al. [29] proposed a similar integrated model based on discrete choice models.

Different to the common features used in the FA problem, Pita et al. [30] considered airports with the slot constrained (SC) in their study. The slot can be defined as the row/time allocation to an aircraft. The aim is to reduce large delays by controlling air traffic. Kenan et al. [31] assigned fleets to flights which have stochastic demand and stochastic fare (SF) based on fare classes. Sherali et al. [32,33] emphasized the importance of considering fuel consumption (FUEL) in airline optimization decisions. Gürkan et al. [34] integrated the cruise time controllability (CTC) of the aircraft using cruise speed adjustments into FA and AR problems. Jamili [35] integrated all three decisions into a single model, but also developed flight programs from the ground up and addressed the issue of fleet type route delays by introducing buffers. Safak et al. [36] took into account passenger connection times, fuel consumption, and CO₂ emission costs associated with cruise speed settings by integrating the three decisions.

Kenan et al. [6] integrated CA into the FS, FA, and AR problems. It highlights that such agreements can significantly affect an airline’s profit by using fewer aircraft, fewer delays, and more flights. They also added optional flight legs (OF) to the model, which can be canceled according to demand changes. Şafak et al. [37] developed a three-stage stochastic programming model for the FS, FA, and AR problems. Cacchiani and Gonzalez [38] aggregated the FS, FA, AR, and CS problems. Xu et al. [39] proposed an integrated model for the FS, FA, and AR problems that allows for spreading delays, on-demand flights, and passenger spills. The spill passenger (SPP) is the number of passengers who cannot be served for any flight.

Kenan et al. [6] and Şafak et al. [36,37,40] have made significant contributions to solving the FA problem. However, Kenan’s study did not detail the costs related to passengers in the model, only considering the delay costs (DC) in the objective function. Şafak et al. [36,37] have discussed costs in detail, and CTC has been used to compensate for the uncertainty of non-cruise time (NCTU). However, different costs have been incurred by increasing the cruise time for new flights or leasing alternatives.

The objective of this study is to maximize the profitability of airline companies, which hinges on the revenue generated from ticket sales to passengers. To efficiently serve a greater number of passengers on each flight, it is imperative to assign appropriate fleet types. Assigning aircraft with a lower capacity than the demand may result in unserved passengers, constituting a spill passenger cost. Furthermore, ensuring timely departures and arrivals is crucial to enhancing passenger satisfaction levels. Delays caused by aircraft idling at airports can lead to disruptions in consecutive flights. These delays, in turn, may cause passengers to miss their connections, especially when the required connection times cannot be met on connecting flights. The impact of these factors on the service quality translates into costs for the companies. Operational costs are incurred by companies for the use of aircraft on each flight, necessitating revenue generation above these costs to ensure profitability. For flights that may not be profitable on their own, companies can opt for less expensive CA to mitigate costs. To the best of our knowledge, there is no study in the literature that handles stochastic demand and stochastic non-cruise times on delays while also considering the impact of CA on the FA problem. This study aims to fill this gap by providing the following main contributions:

• The development of a new third-stage stochastic non-linear programming model that combines FA with operational FS decisions and AR problem using codesharing agreements.
• The consideration of stochastic demand and stochastic non-cruise times simultaneously in the model.
• The proposal of simulation-based metaheuristic (simheuristic) methods for solving the proposed model. The Monte Carlo method is used for simulation, and the Simulated Annealing (SA) and Cuckoo Search (CS) algorithms are used for metaheuristics. Two new solution methods are proposed by integrating simulation and metaheuristics.

  i.

  SA + MC: Simulated Annealing + Monte Carlo

  ii.

  CS + MC: Cuckoo Search + Monte Carlo

Table 1. Features added to the FA.

	Maintaince Planning	Flight Scheduling	Aircraft Routing	Time Window	Stochastic Demand	Spill Passenger Costs	Passenger Connection Time	Delay Costs	Cruise Time Controbility	Stochastic Non-Cruise Time	Codeshare Aggrements
Authors	MP	FS	AR	TW	SD	SPP	PCT	DC	CTC	NCTU	CA
Levin [13]				√
Berge and Hopperstad [11]					√	√
Barnhart et al. [41], Sosnowska and Rolim [42]	√		√
Rexing et al. [14]				√
Yan and Tseng [43]		√
Ahuja et al. [44]	√			√
Lohatepanont and Barnhart [2]		√				√
Li and Wang [45]	√		√
Yan et al. [46]		√
Yan and Chen [47]			√								√
Pilla et al. [26], Jacobs et al. [22]					√
Sherali and Zhu [25], Naumann et al. [27]					√
Dumas et al. [23]					√	√
Haouari et al. [48]	√		√
Sherali et al. [49]	√	√
Zeghal et al. [16]	√		√	√
Sherali et al. [33]		√		√		√
Liang and Chaovalitwongse [50]	√		√
Pita et al. [30]		√		√		√	√	√
Cadarso and Marin [19]		√					√
Atasoy et al. [29]		√				√
Shao et al. [51]	√		√
Gürkan et al. [34]	√	√	√						√
Dong et al. [52]		√
Liu et al. [53]	√		√
Cacchiani and Salazar-González [54]	√		√
Cadarso and Celis [28]		√			√		√
Jamili et al. [35]	√		√				√
Şafak et al. [36]				√					√
Kenan et al. (a) [6]		√	√		√			√			√
Kenan et al. (b) [31]		√			√
Şafak et al. [37]	√	√	√	√	√		√	√	√	√
Şafak et al. [38]	√	√	√	√	√		√	√	√
Wei et al. [55]	√	√	√	√
Cacchiani and Salazar-González [38]	√	√	√	√
Xu et al. [39]	√	√	√			√	√	√
Ahmed et al. [56]	√		√			√
This study	√	√	√	√	√	√	√	√	√	√	√

√: Access.

Table 1. Features added to the FA.

	Maintaince Planning	Flight Scheduling	Aircraft Routing	Time Window	Stochastic Demand	Spill Passenger Costs	Passenger Connection Time	Delay Costs	Cruise Time Controbility	Stochastic Non-Cruise Time	Codeshare Aggrements
Authors	MP	FS	AR	TW	SD	SPP	PCT	DC	CTC	NCTU	CA
Levin [13]				√
Berge and Hopperstad [11]					√	√
Barnhart et al. [41], Sosnowska and Rolim [42]	√		√
Rexing et al. [14]				√
Yan and Tseng [43]		√
Ahuja et al. [44]	√			√
Lohatepanont and Barnhart [2]		√				√
Li and Wang [45]	√		√
Yan et al. [46]		√
Yan and Chen [47]			√								√
Pilla et al. [26], Jacobs et al. [22]					√
Sherali and Zhu [25], Naumann et al. [27]					√
Dumas et al. [23]					√	√
Haouari et al. [48]	√		√
Sherali et al. [49]	√	√
Zeghal et al. [16]	√		√	√
Sherali et al. [33]		√		√		√
Liang and Chaovalitwongse [50]	√		√
Pita et al. [30]		√		√		√	√	√
Cadarso and Marin [19]		√					√
Atasoy et al. [29]		√				√
Shao et al. [51]	√		√
Gürkan et al. [34]	√	√	√						√
Dong et al. [52]		√
Liu et al. [53]	√		√
Cacchiani and Salazar-González [54]	√		√
Cadarso and Celis [28]		√			√		√
Jamili et al. [35]	√		√				√
Şafak et al. [36]				√					√
Kenan et al. (a) [6]		√	√		√			√			√
Kenan et al. (b) [31]		√			√
Şafak et al. [37]	√	√	√	√	√		√	√	√	√
Şafak et al. [38]	√	√	√	√	√		√	√	√
Wei et al. [55]	√	√	√	√
Cacchiani and Salazar-González [38]	√	√	√	√
Xu et al. [39]	√	√	√			√	√	√
Ahmed et al. [56]	√		√			√
This study	√	√	√	√	√	√	√	√	√	√	√

√: Access.

3. Problem Description and Formulation

3.1. Terminology

• Flight leg: Describes a flight of an aircraft from the departure airport to the destination airport.
• Path (itinerary): A sequence of one or more flight legs between a specific origin and destination.
• Fleet type (aircraft type): A certain model of aircraft. All the aircrafts of the same type have the same cockpit configuration, crew qualification requirements, maintenance requirements, and capacity.
• Fleet family (aircraft family): A set of aircraft types, each having the same cockpit configuration and crew qualification requirements. Thus, the same crew can fly any fleet type of the same family.
• Fare class (FC): Available seats on aircrafts are divided into classes according to their fares. The cost of seats in each class is the same. Seat capacities in these classes may be fixed or variable.
• Turnaround time: It is the time required for a fleet type to prepare for flight. This period includes the cleaning of the aircraft, passenger, and baggage movements.
• Cruise time is the time of flight that falls between climb and descent.
• Non-cruise time: It is the sum of the taxi in and out times between an aircraft’s landing and take-off. In other words, it is the time that the aircraft moves on the ground.
• Tax in: It is the duration of the aircraft moving towards the apron after landing.
• Taxi out It is the duration of the aircraft moving from the apron to the take-off field.

3.2. Problem Description

In order to effectively tackle the complex FA problem, it is imperative to consider the interrelated decisions pertaining to FS and AR while also accounting for codesharing agreements. To address the uncertainty associated with demand and non-cruise time, this study presents a novel and integrated approach that models the problem as a three-stage stochastic mixed-integer non-linear programming (TSS_MINLP) framework.

In line with union regulations, the airline industry necessitates operational decisions to be made 2–3 months prior to the scheduled flight. Similarly, codesharing agreements are strategic decisions that require meticulous planning well in advance. However, the uncertainty surrounding the demand for flights during this period poses a significant challenge (as highlighted by Kenan et al. [6]). Thus, this approach first focuses on making decisions regarding FA, routing, and CA for the previously scheduled flights, followed by determining the number of passengers to be accommodated as the flight date approaches, based on these decisions. Consequently, the fleet type route schedule is adjusted accordingly. Thus, the proposed three-stage decision-making model, illustrated in Figure 1, integrates the FA, FS, and AR decisions in a cohesive framework that effectively addresses the uncertainty associated with demand and non-cruise time.

The problem at hand can be divided into three main steps. The first step involves determining the flights that will be operated by the airline and those that will be operated under CA. Based on this, an appropriate fleet type is assigned to each flight, and the routes are created while accounting for the maintenance planning required for the fleet type. Once consecutive flights are completed, the starting and final stations for each aircraft are determined.

Moving on to the second step, the estimated departure time for the airline’s flights is selected from a predetermined time window. This ensures that the departure time is optimized, keeping in mind the operational constraints and other factors that may influence the timing of the flight.

In the final step, the timetable of the routes assigned to each aircraft is updated, and the number of passengers to be carried by the aircraft is determined. This step is critical in ensuring that the aircraft is utilized optimally while adhering to safety standards and operational efficiency. By determining the number of passengers to be carried, the airline can ensure that the aircraft is optimally utilized, and that the entire operation runs smoothly.

Overall, this three-step approach, which involves determining flights and fleet types, optimizing departure times, and updating the timetable while considering passenger count, enables the airline to operate efficiently and effectively.

The problem has several assumptions which are given as follows:

• The analysis of data obtained from the BTS [57] revealed that non-cruise times follows a normal distribution.
• It is assumed that demand follows a uniform distribution with lower and upper parameters covering the minimum and maximum fleet type capacities.
• The departure and arrival times of flights are also pre-determined.
• A time window is used for departure times, and delays that exceed this time window result in additional costs that are included in the objective function.
• Maintenance planning is conducted at either the first or the last airport.
• Codeshare flights are not included in the airline’s routes, and a specific budget is allocated for CA.
• The cruising speed of the fleet type is allowed to vary at certain rates.
• If there are connecting flights for passengers, a minimum time is given for passengers to switch to the next flight, and if this time is not sufficient, passengers miss their connections.

3.3. Mathematical Model

TSS_MINLP model is below:

Sets
S	:	Scenarios  $\mathrm{s}\in \mathrm{S}$
K	:	Fleet types  $\mathrm{k}\in \mathrm{K}$
F	:	Flights  $\mathrm{i}\in \mathrm{F}$
${\mathrm{Y}\mathrm{B}}_{\mathrm{i}}$	:	Flights which have passenger connections from flight  $\mathrm{i}\in \mathrm{F}$
BU	:	Connected flights sets  $\mathrm{i},\mathrm{j}\in \mathrm{F}$
${\mathrm{O}}_{\mathrm{i}}$	:	Flights before flight  $\mathrm{i}\in \mathrm{F}$
${\mathrm{H}}_{\mathrm{i}}$	:	Flights after flight  $\mathrm{i}\in \mathrm{F}$
${\mathrm{S}\mathrm{U}}_{\mathrm{k}}$	:	Last flights of fleet type  $\mathrm{k}\in \mathrm{K}$
${\mathrm{I}\mathrm{U}}_{\mathrm{k}}$	:	First flights of fleet type  $\mathrm{k}\in \mathrm{K}$
${\rm Y}$	:	Different allied airlines   $\mathsf{\upsilon }\in {\rm Y}$
P	:	Codeshare agreements  $\mathrm{p}\in \mathrm{P}$

Parameters
$\left[{\mathrm{T}}_{\mathrm{i}}^{\mathrm{m}\mathrm{i}\mathrm{n}},{\mathrm{T}}_{\mathrm{i}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\right]$	:	Earliest and latest departure time of flight i  $\mathrm{i}\in \mathrm{F}$
${\mathrm{p}\mathrm{r}}_{\mathrm{s}}$	:	Probability of scenario s  $\mathrm{s}\in \mathrm{S}$
${\mathrm{A}\mathrm{N}}_{\mathrm{k}}$	:	Number of current aircraft of fleet type k  $\mathrm{k}\in \mathrm{K}$
${\mathrm{K}\mathrm{a}\mathrm{p}}_{\mathrm{k}}$	:	Number of seats in fleet type k  $\mathrm{k}\in \mathrm{K}$
${\mathrm{D}\mathrm{Y}}_{\mathrm{i}}$	:	Opportunity cost of spilled passengers of flight i  $\mathrm{i}\in \mathrm{F}$
${\mathrm{B}\mathrm{O}\mathrm{S}}_{\mathrm{k}}$	:	Unit idle time cost of fleet type k (per minute)  $\mathrm{k}\in \mathrm{K}$
${\mathrm{K}\mathrm{Y}}_{\mathrm{i}}$	:	Cost per passenger for miss-connected passengers on flight i  $\mathrm{i}\in \mathrm{F}$
${\mathrm{D}\mathrm{C}}_{\mathrm{i}}$	:	Per minute delay cost of flight i  $\mathrm{i}\in \mathrm{F}$
${\mathrm{f}\mathrm{a}\mathrm{r}\mathrm{e}}_{\mathrm{i}}$	:	Flight revenue per passenger  $\mathrm{i}\in \mathrm{F}$
${\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}}_{\mathrm{i}}$	:	Flight operating cost per passenger  $\mathrm{i}\in \mathrm{F}$
${\mathsf{\mu }}_{\mathrm{i}\mathrm{s}}$	:	Demand of flight i in scenario s  $\mathrm{s}\in \mathrm{S}, \mathrm{i}\in \mathrm{F}$, stochastic
${\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}$	:	Non-cruise time of flight i in scenario s  $\mathrm{s}\in \mathrm{S}, \mathrm{i}\in \mathrm{F}$, stochastic
$\left[{\mathrm{C}\mathrm{T}}_{\mathrm{i}\mathrm{k}}^{\mathrm{l}},{\mathrm{C}\mathrm{T}}_{\mathrm{i}\mathrm{k}}^{\mathrm{u}}\right]$	:	Time window for cruise time of flight i of fleet type k  $\mathrm{i}\in \mathrm{F}, \mathrm{k}\in \mathrm{K}$
${\mathrm{T}\mathrm{A}}_{\mathrm{i}\mathrm{k}}$	:	Turnaround time required to prepare fleet type k after flight i $\mathrm{i}\in \mathrm{F}, \mathrm{k}\in \mathrm{K}$
${\mathrm{C}\mathrm{P}}_{\mathrm{i}\mathrm{j}}$	:	Transit time for connected passengers between flights i, j  $\mathrm{i}\in \mathrm{J}, \mathrm{j}\in {\mathrm{Y}\mathrm{B}}_{\mathrm{i}}$
${\mathrm{p}\mathrm{a}\mathrm{s}\mathrm{s}}_{\mathrm{i}\mathrm{j}}$	:	Number of passengers from flight i connected to flight j  $\mathrm{i}\in \mathrm{J}, \mathrm{j}\in {\mathrm{Y}\mathrm{B}}_{\mathrm{i}}$
${\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{e}}_{\mathsf{\upsilon }\mathrm{p}\mathrm{i}}$	:	The cost per passenger in contract type p  $\mathrm{p}\in \mathrm{P}, \mathsf{\upsilon }\in {\rm Y}, \mathrm{i}\in \mathrm{F}$
${\mathrm{C}\mathrm{a}\mathrm{p}}_{\mathsf{\upsilon }\mathrm{p}\mathrm{i}}$	:	Capacity provided under contract type p on flight leg i, $\mathrm{p}\in \mathrm{P}, \mathsf{\upsilon }\in {\rm Y}, \mathrm{i}\in \mathrm{F}$
BUD	:	Available budget for codeshare agreements
θ	:	The ratio of total codeshare flight capacities to the total capacity of the airline
M	:	A big number

Decision Variables
${\mathrm{x}}_{\mathrm{i}\mathrm{j}\mathrm{k}}^1$	:	$\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\mathrm{If}\mathrm{flights}\mathrm{j}\mathrm{comes}\mathrm{after}\mathrm{flight}\mathrm{i}\mathrm{and}\mathrm{both}\mathrm{are}\mathrm{performed}\mathrm{with}\mathrm{fleet}\mathrm{type}\mathrm{k}\hspace{1em}\hspace{1em}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{j}\in \mathrm{F}, \mathrm{i} \in {\mathrm{O}}_{\mathrm{j}}, k\in K\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$
${\mathrm{f}}_{\mathrm{i}\mathrm{k}}^1$	:	$\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\mathrm{If}\mathrm{flight}\mathrm{i}\mathrm{is}\mathrm{the}\mathrm{first}\mathrm{flight}\mathrm{with}\mathrm{fleet}\mathrm{type}\mathrm{k}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{i}\in \mathrm{F}, \mathrm{k}\in K\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$
${\mathrm{e}}_{\mathrm{i}\mathrm{k}}^1$	:	$\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\mathrm{If}\mathrm{flight}\mathrm{i}\mathrm{is}\mathrm{the}\mathrm{last}\mathrm{flight}\mathrm{with}\mathrm{fleet}\mathrm{type}\mathrm{k}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{i}\in \mathrm{F}, \mathrm{k}\in \mathrm{K}\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$
${\mathrm{q}}_{\mathsf{\upsilon }\mathrm{p}\mathrm{i}}^1$	:	$\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\mathrm{If}\mathrm{the}\mathrm{flight}\mathrm{leg}\mathrm{i}\mathrm{is}\mathrm{performed}\mathrm{with}\mathrm{the}\mathrm{aircraft}\mathrm{in}\mathrm{the}\mathrm{contract}\mathrm{p}\mathrm{of}\mathrm{company} \mathsf{\upsilon }\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\upsilon }\in {\rm Y}, \mathrm{i} \in \mathrm{F}, \mathrm{p} \in \mathrm{P}\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$
${\mathrm{a}}_{\mathrm{i}}^2$	:	Announced departure time of flight i  $\mathrm{i}\in \mathrm{F}$
${\mathrm{b}}_{\mathrm{i}\mathrm{s}}^3$	:	Actual departure time of flight i in scenario s  $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}$
${\mathrm{c}}_{\mathrm{i}\mathrm{s}}^3$	:	Actual arrival time of flight i in scenario s  $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}$
${\mathrm{d}}_{\mathrm{i}\mathrm{k}\mathrm{s}}^3$	:	Cruise time of flight i with fleet type k in scenario s  $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}, \mathrm{k}\in \mathrm{K}$
${\mathrm{I}\mathrm{T}}_{\mathrm{i}\mathrm{k}\mathrm{s}}^3$	:	Idle time of flight i with fleet type k in scenario s  $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}, \mathrm{k}\in \mathrm{K}$
${\mathrm{d}\mathrm{e}\mathrm{l}}_{\mathrm{i}\mathrm{s}}^3$	:	Delay time of flight i in scenario s  $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}$
${\mathsf{\pi }}_{\mathrm{i}\mathrm{s}}^3$	:	Number of accepted passengers in flight leg i under scenario s $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}$
${\mathsf{\pi }\mathrm{c}}_{\mathrm{i}\mathrm{s}}^3$	:	Number of accepted passengers in flight leg i under scenario s in codeshare $\mathrm{i}\in \mathrm{F}, \mathrm{s}\in \mathrm{S}$
${\mathrm{w}}_{\mathrm{i}\mathrm{j}\mathrm{s}}^3$	:	$\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\mathrm{In}\mathrm{scenario}\mathrm{s}, \mathrm{if}\mathrm{passengers}\mathrm{on}\mathrm{flight}\mathrm{i}\mathrm{miss}\mathrm{flight}\mathrm{j}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{i}\in \mathrm{J}, \mathrm{j} \in {YB}_{\mathrm{i}}, \mathrm{s} \in \mathrm{S},\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$

Maximize

$$={\displaystyle \sum _{s\in S}}{pr}_s\left[\begin{array}{c}\left({\displaystyle \sum _{i\in F}}\left({\pi }_{is}^3+{\pi c}_{is}^3\right)\ast {fare}_i\right)-\left({\displaystyle \sum _{i\in F}}{\pi }_{is}^3\ast {cost}_i\right)\\ -\left({\displaystyle \sum _{iϵF}}{\displaystyle \sum _{\mathsf{\upsilon }\in {\rm Y}}}{\displaystyle \sum _{pϵP}}{code}_{\mathsf{\upsilon }pi}\ast {Cap}_{\mathsf{\upsilon }pi}\ast q_{\mathsf{\upsilon }pi}^1\right)\\ -\left({\displaystyle \sum _{iϵF}}{\displaystyle \sum _{\mathsf{\upsilon }\in {\rm Y}}}{\displaystyle \sum _{pϵP}}\max (0,{\mu }_{is}-{Cap}_{\mathsf{\upsilon }pi})\ast q_{\mathsf{\upsilon }pi}^1\ast \left({fare}_i-{code}_{\mathsf{\upsilon }pi}\right)\right)\\ -\left({\displaystyle \sum _{iϵF}}{\displaystyle \sum _{kϵK}}\max \left(0,{\mu }_{is}-{Kap}_k\right)\ast \left({\displaystyle \sum _{jϵO_i}}{x}_{jik}^1+f_{ik}^1\right)\ast {DY}_i\right)\\ -\left({\displaystyle \sum _{iϵF}}{\displaystyle \sum _{kϵK}}{IT}_{iks}^3\ast {BOS}_k\right)-\left({\displaystyle \sum _{iϵF}}{DC}_i\ast {del}_{is}^3\right)-\left({\displaystyle \sum _{iϵF}}{\displaystyle \sum _{jϵ{YB}_i}}{KY}_j\ast {pass}_{ij}\ast w_{ijs}^3\right)\end{array}\right]$$

(1)

Subject to

$$\sum _{k\in K}\left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right)+\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{p\in P}{q}_{\mathsf{\upsilon }pi}^1=1 \forall i\in F$$

(2)

$$\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1-\sum _{jϵH_i}{x}_{ijk}^1-e_{ik}^1=0 \forall i\in F,\forall k\in K$$

(3)

$$\sum _{i\in F}{f}_{ik}^1\le {AN}_k \forall k\in K$$

(4)

$$f_{ik}^1=0 \forall i\in F\backslash {IU}_k,\forall k\in K$$

(5)

$$e_{ik}^1=0 \forall i\in F\backslash {SU}_k,\forall k\in K$$

(6)

$$\sum _{iϵF}\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{pϵP}{code}_{\mathsf{\upsilon }pi}\ast {Cap}_{\mathsf{\upsilon }pi}\ast q_{\mathsf{\upsilon }pi}^1\le BUD$$

(7)

$$\sum _{iϵF}\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{pϵP}{Cap}_{\mathsf{\upsilon }pi}\ast q_{\mathsf{\upsilon }pi}^1\le \theta \ast \sum _{k\in K}{AN}_k\ast {Kap}_k$$

(8)

$$\sum _{iϵF}\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{pϵP}{q}_{\mathsf{\upsilon }pi}^1\le \sum _{iϵF}\sum _{kϵK}\left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right)$$

(9)

$$\left[1-\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{p\in P}{q}_{\mathsf{\upsilon }pi}^1\right]\ast T_i^{min}\le {\mathrm{a}}_{\mathrm{i}}^2\le \left[1-\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{p\in P}{q}_{\mathsf{\upsilon }pi}^1\right]\ast T_i^{max} \forall i\in F$$

(10)

$$-{\mathrm{a}}_{\mathrm{i}}^2+b_{is}^3\ge 0 \forall i\in F, \forall s\in S$$

(11)

$$b_{is}^3-{\mathrm{a}}_{\mathrm{i}}^2-{del}_{is}^3\le 0 \forall i\in F, \forall s\in S$$

(12)

$$b_{is}^3+\sum _{k\in K}\left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right)\ast {\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}+\sum _{k\in K}{d}_{iks}^3=c_{is}^3 \forall i\in F,\forall s\in S$$

(13)

$$If x_{ijk}^1=1, b_{js}^3-b_{is}^3={\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}+{\mathrm{T}\mathrm{A}}_{\mathrm{i}\mathrm{k}}+d_{iks}^3+{IT}_{iks}^3 \forall \left(i,j\right)\in BU, \forall s\in S,\forall k\in K$$

(14)

$$-{CT}_{ik}^l\left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right)+d_{iks}^3\ge 0 \forall i\in F, \forall s\in S,\forall k\in K$$

(15)

$${CT}_{ik}^u\left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right)-d_{iks}^3\ge 0 \forall i\in F, \forall s\in S,\forall k\in K$$

(16)

$$\begin{gathered} b_{js}^3-b_{is}^3-\sum _{k\in K}{d}_{iks}^3+M\ast w_{ijs}^3\ge \left[1-\sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{p\in P}{q}_{\mathsf{\upsilon }pi}^1\right]\ast \left({\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}+{CP}_{ij}\right) \\ \forall i\in F , \forall j\in {YB}_i ,\forall s\in S \end{gathered}$$

(17)

$${\pi }_{is}^3\le \sum _{k\in K}{Kap}_k\ast \left(\sum _{jϵO_i}{x}_{jik}^1+f_{ik}^1\right) \forall i\in F,\forall s\in S$$

(18)

$${\pi c}_{is}^3\le \sum _{\mathsf{\upsilon }\in {\rm Y}}\sum _{pϵP}{Cap}_{\mathsf{\upsilon }pi}\ast q_{\mathsf{\upsilon }pi}^1 \forall i\in F, \forall s\in S$$

(19)

$${\pi }_{is}^3\le {\mu }_{is} \forall i\in F,\forall s\in S$$

(20)

$${\pi c}_{is}^3\le {\mu }_{is} \forall i\in F, \forall s\in S$$

(21)

$$x_{ijk}^1 \in \left\{\mathrm{0,1}\right\} \forall i\in F,\forall j\in O_i, \forall k\in K$$

(22)

$$f_{ik}^1 \in \left\{\mathrm{0,1}\right\} \forall i\in F, \forall k\in K$$

$$e_{ik}^1 \in \left\{\mathrm{0,1}\right\} \forall i\in F, \forall k\in K$$

$$q_{\mathsf{\upsilon }pi}^1 \in \left\{\mathrm{0,1}\right\} \forall i\in F, \forall \mathsf{\upsilon }\in {\rm Y},\forall p\in P$$

$$\begin{gathered} {\mathrm{a}}_{\mathrm{i}}^2\ge 0 \forall i\in F \\ {\pi }_{is}^3\ge 0 \forall i\in F , \forall s\in S \\ {\pi c}_{is}^3\ge 0 \forall i\in F , \forall s\in S \\ {IT}_{iks}^3\ge 0 \forall i\in F , \forall k\in K,\forall s\in S \end{gathered}$$

$$\begin{gathered} {del}_{is}^3\ge 0 \forall i\in F, \forall s\in S \\ {w}_{ijs}^3\in \left\{\mathrm{0,1}\right\} \forall i\in F,\forall j\in {YB}_i,\forall s\in S \end{gathered}$$

Equation (1) aims to maximize the expected profit of an airline, which is calculated as the difference between the income from passengers and the associated costs. The costs include flight cost per passenger, codeshare costs, cost of spilled passengers, fleet type idle time cost, delay cost, and passenger misconnection costs. Equations (2)–(22) are constraints. The decisions of the first stage are determined by Equations (2)–(9). Equation (2) ensures that each planned flight leg is covered by exactly one fleet type, while Equation (3) maintains the route structure of the network. Equation (4) specifies that the total number of aircraft which is used should not exceed the number of aircraft available. Equations (5) and (6) identify the first and last airports of each fleet type. These airports serve as maintenance stations for the aircraft. Equation (7) ensures that the total investment in CA between operators and flight legs does not exceed the allocated budget. The airline being considered is a marketing carrier that operates its own flights and may also serve as the operator for other flights. Hence, Equation (8) specifies that the ratio of the total number of seats purchased to the total capacity of all available aircraft should not exceed the value of (θ). Equation (9) specifies that the number of codeshare flights should not exceed the number of the company’s own flights. The second-stage decision is determined by Equation (10), which ensures that the announced departure time of the airline’s own flights falls between determined lower and upper limits. Equations (11)–(21) determine the third-stage decisions of the model. Equations (11) and (12) allow for delays between the announced departure times of the first stage and the actual departure times. Since these delays can impact other flights on the route, they are considered as costs in the objective function. Equation (13) calculates the arrival times for each flight based on the actual departure times. Equation (14) determines the departure times of successive flights on a route, as well as the minimum idle times of fleet type. Equations (15) and (16) impose upper and lower cruise time limits for each flight. Equation (17) accounts for passengers who miss their connecting flights due to the absence of sufficient transit time, which creates a cost for the company. Equations (18) and (19) specify that the number of passengers admitted to a flight leg must not exceed the capacity of the aircraft or the capacity in the selected contract. Similarly, Equations (20) and (21) ensure that the number of passengers accepted on a flight leg does not exceed the demand for the flight leg. Finally, Equation (22) defines the domain of all variables.

The proposed model is subjected to a range of diverse scenarios, characterized by distinct levels of passenger demand and non-cruise time of flights. The resulting third-stage decisions and corresponding objective function values exhibit significant variability across these scenarios. To effectively account for these real-world complexities, the model employs a probabilistic weighting approach. Specifically, the outcomes for each scenario are multiplied by their respective probabilities, and subsequently summed to derive the expected value of the objective function. This approach enables the model to incorporate and adapt to the diverse conditions encountered in practical settings.

4. Solution Methodology

Optimization packages are limited in their ability to solve the TSS_MINLP, especially when the number of flights and scenarios is increased. Therefore, the development of heuristic or metaheuristic methods is necessary [34,35,36,38,39,55].

Simheuristics are an extension of metaheuristics for solving optimization problems under uncertainty. They have been widely applied to different problems such as production planning, supply chain design, flexible flow shop scheduling, resource allocation, transportation, and healthcare scheduling [58,59,60,61,62,63]. Simheuristics combine metaheuristics with simulation methods, allowing for stochastic variables to be considered in mathematical models [64]. This means that a simheuristic algorithm includes a special simulation for stochastic variables, making them efficient and effective in overcoming stochastic optimization problems [65]. In the proposed solution method, scenarios with varying demand and non-cruise times are generated using the Monte Carlo simulation method.

4.1. Simulated Annealing

Simulated annealing (SA) is a search algorithm that is inspired by the physical annealing process of a solid material. Physical annealing involves gradually cooling a heated solid to its minimum energy state [66]. SA begins by performing a broad search of the solution space, accepting even suboptimal solutions with a high probability. As the algorithm progresses, it performs a more focused search around the current solution, with the goal of finding the global optimum. At each iteration, SA compares the current solution to a neighboring solution. If the neighboring solution is better, it is accepted as the new current solution. If the neighboring solution is worse, the algorithm may still accept it with a certain probability, to avoid getting stuck at a local minimum. The best available solution is also retained. The acceptance probability is based on the change in profit and the temperature parameters. This process is repeated at each temperature level, with the temperature gradually decreasing according to a predetermined cooling schedule [67].

SA is a stochastic search algorithm used to reach the global optimum hidden among local optima. It has the great advantage of having a diverse and deep neighborhood structure. It is also a simple and fast algorithm with few parameters [68]. In the literature, strong results have been obtained in most applications [69]. Therefore, this study uses SA as a metaheuristic method, integrated with the Monte Carlo simulation (MC). The flow diagram of the simheuristic method (SA + MC) is shown in Figure 2, and the detailed pseudocode is given in Algorithm 1. The flowchart is divided into smaller rectangles, numbered, and drawn with dashed lines to aid in explanation. The integration of SA and MC allows for the efficient and effective optimization of complex problems with stochastic variables.

Figure 2 provides an overview of the proposed methodology, which comprises several iterative stages. In the first red rectangle, an initial solution is generated by creating routes that are suitable for flights. Next, the fleet type with the closest capacity to the flight with the highest expected demand on each route is assigned to that route. The resulting assignment and routing decisions are evaluated using a short Monte Carlo simulation (e.g., 100 runs) to calculate the expected objective function value of the problem, as shown in the second red rectangle.

Subsequently, the SA algorithm is initiated with the initial solution obtained in the first red rectangle. At each iteration, a neighboring solution is generated using the algorithm presented in Algorithm 2, as shown in the third red rectangle. The neighboring solution assigns different fleet types to some flights and creates new routes based on this assignment. The expected objective function value is then calculated in a short Monte Carlo simulation of the newly generated neighboring solution, as shown in the fourth red rectangle.

In the fifth red rectangle, the SA algorithm updates the best solution and best objective function value based on the newly generated neighboring solution. Once the stopping criteria are met, the expected objective function value of the best solution in the long Monte Carlo simulation (e.g., 1000 runs) is updated, as shown in the sixth red rectangle.

Algorithm 1: SA + MC Algorithm

Input flight informations Create initial solution $s_0 \in S:$   Create feasible routes for flights Generate random values $\left({\mathsf{\mu }}_{\mathrm{i}\mathrm{s}}, {\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}\right)$ from appropriate distributions for scenarios   Calculate the expected demand value for each flight Assign capacity airplanes that match these demand values Calculate the objective function value $f\left(s_0\right)$ with short Monte Carlo simulation of the initial solution Set an initial temperature: $T>0;$                  $s=s_0; f\left(s\right)=f\left(s_0\right);$                 $s_{iyi}=s_0; f\left(s_{iyi}\right)=f\left(s_0\right);$ Repeat                    $\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}=1$ Repeat Generate a neighbor solution $s^{\prime }$ of S Create routes based on new assignment Calculate the objective function value $f\left(s^{\prime }\right)$ with a short Monte Carlo simulation of $s^{\prime }$ solution $\Delta =f\left(s^{\prime }\right)-f\left(s\right);$ If $\Delta \ge 0$ then $s=s^{\prime };$ If $\Delta <0$ then generate a random number u in the range (0,1) If $u<e^{\frac{\Delta }{T}}$ then $s=s^{\prime };$ If $f\left(s^{\prime }\right)<f\left(s_{iyi}\right)$ then $s_{iyi}=s^{\prime }$                   $\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}=\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}+1;$ Until $count>M$ (number of neighbors to search)                     $i=i+1;$                   $T\left(i\right)=\alpha \ast T\left(i-1\right);$ Until (until the stopping criteria is met) $s_{iyi}$ best solution Calculate objective function value f(s_iyi) with long Monte Carlo simulation

Algorithm 2: Neighbor Solution Search Algorithm

$\mathrm{f}\mathrm{l}\mathrm{t}=0,$ y, generate a random value between 0–1 If $y\le 0.3$ then  $\mathrm{For} (i=0; i<3; i++)$ Assign fleet type assigned to randomly selected flight in the same solution to flight $(flt+i)$ in solution End for If not    For $(i=0; i<3; i++)$    Assign random fleet type to flight $(flt+i)$ in solution    End for End if $flt=flt+1$

4.2. Cuckoo Search Algorithm (CS)

The CS algorithm draws inspiration from the aggressive breeding techniques employed by cuckoos. Cuckoos lay their eggs in the nests of other birds that are randomly selected. If the host bird discovers the cuckoo egg, it may either remove the egg or abandon the nest entirely. The eggs that survive form the next generation [70].

The study by Rakesh and Mahesh [71] describes the successful application of the nature-inspired cuckoo algorithm, with various variations, to comprehensive global optimization problems. CS has been shown to outperform other metaheuristics in terms of speed and accuracy in solving the traveling salesman problem (TSP) [72]. Additionally, it has demonstrated a superior ability to find higher quality solutions for constrained optimization problems than genetic algorithm and particle swarm optimization. Its stochastic nature makes it particularly effective in exploring the search space, allowing for its successful application to highly constrained optimization problems such as job scheduling, multi-objective scheduling, flow shop scheduling, and TSP. However, in such problems, some modifications may be required to produce new solutions with Levy flights [73,74,75]. New solutions can also be produced with the help of heuristics without using the Levy distribution. Gao et al. [76]’s study contains many examples. In this study, improved new solutions with local search heuristics (Algorithm 3) are produced without using levy flights. The local search (Algorithm 4) is also used in the SA + MC algorithm.

The CS algorithm initially starts with a certain number of candidate solutions, then evaluates new solutions and repeats a certain number of bad solutions. The search process in this solution space shows that it has potential to better balance exploitation and exploration. The algorithm has also fewer adjustable parameters [77]. The reason for using CS is that it has these advantages along with its performance in the literature.

The flow diagram of the resulting simheuristic method (CS + MC) is shown in Figure 3. The detailed pseudo code of this flowchart is given in Algorithm 3. The flowchart is divided into small rectangles drawn with dashed lines and numbered to aid in explanation.

Figure 3 illustrates the multi-stage process of the proposed methodology for solving the airline profit maximization problem using the CS algorithm. In the first red rectangle, an initial solution equal in size to the population is generated, following the same approach as is depicted in Figure 2.

In the second red rectangle, the expected objective function value of each initial solution is calculated using a short Monte Carlo simulation (e.g., 100 runs). The initial population is then fed into the CS algorithm, as shown in the third red rectangle.

At each iteration, a randomly selected solution is updated using the local search algorithm outlined in Algorithm 4, as shown in the third red rectangle. The updated solution is then compared to another randomly selected solution from the population, and the best solution and objective function value are updated accordingly. This process is iteratively repeated until the stopping criterion is satisfied.

Finally, the objective function value of the best solution is updated using a long Monte Carlo simulation (e.g., 1000 runs), as shown in the fourth red rectangle. This comprehensive approach ensures that the optimal solution is robust and reliable, accounting for the inherent complexities and uncertainties of real-world airline operations.

Algorithm 3: CS + MC Algorithm

Input flight informations Generate P initial populations $x_i$ (i = 1, …, P)   Create feasible routes for flights Generate random values $\left({\mathsf{\mu }}_{\mathrm{i}\mathrm{s}}, {\mathsf{\eta }}_{\mathrm{i}\mathrm{s}}\right)$ from appropriate distributions for scenarios   Calculate the expected demand value for each flight Assign capacity airplanes that match these demand values Calculate the objective function value $f\left(x_i\right)$ with short Monte carlo simulation of $x_i$ solution While $(t<Maximum iteration)$ do Take a random cuckoo $x_i$ Local search $x_i^{\prime }$ and $f\left(x_i^{\prime }\right)$ Randomly select a nest (let it be $x_j$) in P nests If $(f\left(x_i^{\prime }\right)>f\left(x_j\right))$ then Replace $x_j$ with new solution end if Abandon the worst nests with α ratio of $p_{\alpha }$ and create new ones   Create routes for flights   Assign a random fleet type to each route Calculate the objective function value with a short Monte Carlo simulation of each new solution Keep best solutions (nests with quality solution) Sort the solutions and find the current best solution end while Calculate objective function value with long Monte Carlo simulation of best solution

Algorithm 4: Local Search Algorithm

For $(fl=0; fl<flightnumber-2; fl++)$ y, generate a random value between 0–1 If $y\le 0.3$ then   For $(i=0; i<3; i++)$ Assign fleet type assigned to randomly selected flight in the same solution to flight $(fl+i)$ in solution $x_i$ Create routes based on new assignment Calculate the objective function value $f\left(x_i^{\prime }\right)$ with short Monte Carlo simulation of new solution $\left(x_i^{\prime }\right)$ End for If not     For $(i=0; i<3; i++)$     Assign random fleet type to flight $(fl+i)$ in solution $x_i$ Create routes based on new assignment Calculate the objective function value $f\left(x_i^{\prime }\right)$ with short Monte Carlo simulation of new solution $\left(x_i^{\prime }\right)$    End for End if End for

4.3. Small Example

This section aims to demonstrate the validation and verification of the proposed solution methods through a small sample problem and its corresponding solution. The sample problem includes six flights, two fleet types, two codeshare agreements, and five airports. These airports are Los Angeles (LAX), Dallas Fort Worth (DFW), Atlanta (ATL), Miami (MIA), and Kahului (OGG). The parameters of the flights are given in

Table 2. Flight parameters of small example.

							Min. Departure Time		Max. Departure Time
Flights (i)	Departure Airports	Arrival Airports	Departure Time		Arrival Time		${\mathbf{T}}_{\mathbf{i}}^{\mathbf{m}\mathbf{i}\mathbf{n}}$ (min)		${\mathbf{T}}_{\mathbf{i}}^{\mathbf{m}\mathbf{a}\mathbf{x}}$ (min)
1	LAX	ATL	08:00		12:30		470		490
2	ATL	LAX	13:10		17:50		780		800
3	LAX	DFW	18:40		21:55		1110		1130
4	LAX	MIA	05:00		10:10		290		310
5	MIA	LAX	11:30		16:50		680		700
6	LAX	OGG	18:00		23:40		1070		1090
	Spilled Passengers Cost		Misconnected Passenger Cost		Aircraft Delay Cost		Flight Fare per Passenger		Flight Operating Cost per Passenger
Flights (i)	${\mathrm{D}\mathrm{Y}}_{\mathrm{i}}$ (USD)		${\mathrm{K}\mathrm{Y}}_{\mathrm{i}}$ (USD)		${\mathrm{D}\mathrm{C}}_{\mathrm{i}}$ (USD)		${\mathrm{f}\mathrm{a}\mathrm{r}\mathrm{e}}_{\mathrm{i}}$ (USD)		${\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}}_{\mathrm{i}}$ (USD)
1	102		50		4.85		255		153
2	102		50		7.03		255		153
3	88		50		4.29		220		132
4	163.6		50		3.84		409		245.4
5	124		50		7.21		310		186
6	89.6		50		5.51		224		134
	Codesharing Costs per Passenger		Codesharing Capacity		Min. Cruise Time		Max. Cruise Time		Turnaround Time
	${\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{e}}_{\mathsf{\upsilon }\mathrm{p}\mathrm{i}}$ (USD)		${\mathrm{C}\mathrm{a}\mathrm{p}}_{\mathsf{\upsilon }\mathrm{p}\mathrm{i}}$		${\mathrm{C}\mathrm{T}}_{\mathrm{i}\mathrm{k}}^{\mathrm{l}}$ (min)		${\mathrm{C}\mathrm{T}}_{\mathrm{i}\mathrm{k}}^{\mathrm{u}}$ (min)		${\mathrm{T}\mathrm{A}}_{\mathrm{i}\mathrm{k}}$ (min)
Flights (i)	p = 1	p = 2	p = 1	p = 2	k = 1	k = 2	k = 1	k = 2	k = 1	k = 2
1	76.5	65.03	80	110	213	213	250	250	42	39
2	77	65	80	110	221	221	260	260	46	43
3	66	56.1	80	110	149	149	175	175	44	42
4	122.7	104.3	80	110	247	247	290	290	50	47
5	93	79.05	80	110	255	255	300	300	46	43
6	67.2	57.12	80	110	272	272	320	320	26	25

. There are two fleet types in the problem: B737-800 and A321-200. The parameters for the fleet type are also given in

Table 3. Fleet type parameters.

	Number of Aircrafts	Capacity of Fleet Type	Aircraft Idle Time Costs
Fleet type (k)	${\mathrm{A}\mathrm{N}}_{\mathrm{k}}$	${\mathrm{K}\mathrm{a}\mathrm{p}}_{\mathrm{k}}$	${\mathrm{B}\mathrm{O}\mathrm{S}}_{\mathrm{k}}$ (USD/min)
1	1	172	40
2	1	187	55

.

The problem has been solved for two scenarios by using the GAMS/BARON optimization package and the proposed solution methods. The solutions, presented in

Table 4. Solution results of the sample problem.

	GAMs/BARON	SA + MC	CS + MC
Objective function values (USD)	89,422.62	89,422.62	89,422.62
Assignments	1-codeshare  4-A321-200 2-codeshare  5-A321-200 3-codeshare  6-A321-200
Route	4-5-6
Stations	Codeshare: LAX-ATL Codeshare: ATL-LAX Codeshare: LAX-DFW Route: LAX-MIA-LAX-OGG

, show that the GAMS/BARON and proposed solution methods have obtained the same assignments and routes for each scenario. According to the solution, Flights 1, 2, and 3 are codeshare flights and are operated with aircraft under the second type contract. Therefore, these flights are not included in the airport’s route plan. Flights 4, 5, and 6 are operated by the airline’s second type of aircraft, specifically, the A321-200. The expected objective function value has been obtained by averaging the objective functions of the scenarios. The results of the SA + MC and CS + MC are consistent with GAMS result. Therefore, the validation and verification of the proposed methods have been successfully completed. The detailed plan for each scenario that includes the actual departure time, actual arrival time, cruise time, idle time, delay time, and non-cruise time are listed in

Table 5. Solutions of the scenarios.

	Flights	Actual Departure Time	Actual Arrival Time	Cruise Time (min)	İdle Time (min)	Delay Time (min)	Non-Cruise Time (min)	Demand	Passenger Number
First scenario results	1							147	110
	2							74	74
	3							220	110
	4	05:44	10:58	290	0	34	24	129	129
	5	11:45	17:07	300	0	5	22	201	187
	6	17:50	22:46	272	0		24	172	172
Second scenario results	1							83	83
	2							136	110
	3							65	65
	4	05:45	10:57	290	0	35	22	138	138
	5	11:44	17:07	300	0	4	23	159	159
	6	17:50	22:43	272	0		21	170	170

.

5. Computational Experiments

5.1. Data Descriptions

Four test problems have been generated using data obtained from the “Airline On-Time Performance Data” open access database (BTS) [57]. The flight information has been filtered to include only flights departing and landing at Los Angeles City airport (LAX) from American Airlines. The resulting test problems have 30, 78, 127, and 180 flights, respectively, and include six fleet types.

Table 6. Properties of fleet types.

Fleet Type	B737-800	A321-200	A319-100	A321-NEO	B787-8	ERJ-175
Capacity	172	187	128	196	234	76
Base Turntime	25.47	25.47	23.87	31.05	42.39	16.47
Idle Time Costs (USD/min)	140	142	136	144	147	125

provides general information about the fleet types used in the test problems.

In the context of codeshare flights, the model assumes the existence of two contract types, which entail the operating carrier receiving 30% and 23% of the ticket price, respectively. The passenger capacities associated with these contracts are set at 56 and 80. To obtain the codeshare ticket prices, the model leverages data from the airline’s website. It should be noted that the total passenger capacity on codeshare flights is subject to the constraint that it cannot exceed one third of the airline’s overall aircraft capacity.

It is worth mentioning that the model relies on a number of parameters that have been previously determined by Şafak et al. [36,37]. These parameters are critical to ensuring the validity and reliability of the model’s outputs:

• A 10 min tolerance period is added to the departure times.
• The turnaround time of each fleet type is multiplied by the complexity coefficient of the airport for each flight and fleet type to obtain the turnaround time for that flight.
• The connection time for the passengers is set at 30 min.
• The highest and lowest cruise times are calculated by subtracting 20 min from the total flight time. The lowest cruise times are then compressed by 15% of the highest cruise time.
• The profit from one passenger is defined as the basic spill cost of a lost passenger. The passenger spill cost for each flight has been calculated by multiplying the basic spill cost by the airport complexity coefficient.
• The cost of a missed connection is set at USD 50 per passenger.

In order to solve the problem with 1000 scenarios, 30 flights, and 6 fleet types, the corresponding TSS_MINLP formulation has 704,532 constraints, 518,577 variables, and 3,122,383 non zeroes after the pre-solve is completed. This calculation also describes the complexity of the algorithm. As the size of the problem increases, it becomes more difficult to solve. Therefore, simheuristic solution methods have been used to solve the test problems. Based on experimental studies, the following parameters have been used for the SA algorithm: an initial temperature of 500, a final temperature of 10, a search for 100 neighbors in each step, and a temperature reduction ratio of 0.99. For the CS algorithm, the population size is set to 20, the rate of changing bad solutions is set to 0.2, and the maximum number of iterations is set to 500.

5.2. Computational Results

The proposed methods for solving large-scale problems have been implemented using a C# program and utilized to solve the test problems. Their performances have been compared by using the same stochastic values. Each test problem has been executed with 100 scenarios in a short simulation and 1000 scenarios in a long simulation, with the number of aircraft and flights per fleet type varying based on the problem size. The test problems have been solved using the proposed solution methods for the same total run times, and the best results have been given in

Table 7. Comparison of results of test problems.

				SA + MC			CS + MC
Problems	Scenario Numbers	Flight Numbers	Aircraft Numbers	Obj. Function Values (USD)	Code Flight Numbers	Times (s)	Obj. Function Values (USD)	Code Flight Numbers	Times (s)
1	100	30	6	319,314.3	0	209.26	322,102.7	0	216.01
2	100	78	14	767,161.46	25	665.96	874,041.44	2	758.32
3	100	127	20	1,327,546.52	28	1467.69	1,436,292.3	12	2406.63
4	100	180	25	1,731,738.91	56	2097.55	1,947,504.96	12	3432.55

. The number of algorithm runs during these periods and the results obtained in each run are presented in Figure 4a–d. The results indicate that CS + MC achieved better profits than SA + MC, possibly due to the CS algorithm’s population-based approach, which expands the search space with more diverse solutions. Furthermore, unlike in Şafak’s study [37], where scenarios have been grouped and tested, this study tested many scenarios using metaheuristics and obtained fast and effective results.

5.3. Effect of Codeshare Flights

To analyze the impact of the codesharing strategy, test problem 4 (180 flights) has been solved again without the codesharing option, and the obtained solutions are presented in

Table 8. Results of proposed methods for 180 flights without codeshare flights.

	SA + MC			CS + MC
Run	Profit (USD)	Route Number	Deadhead Flights	Profit (USD)	Route Number	Deadhead Flights
1	1,270,180.8	44	21	1,605,883.39	46	14
2	1,048,609.36	44	21	1,553,091.89	45	15
3	1,138,399.14	44	21	1,572,453.71	46	10
4	1,000,286.11	45	17	1,477,181.86	46	12
5	1,240,760.95	45	17	1,618,113.63	45	17
6	1,282,708.31	44	20	1,479,553.91	45	15
7	1,183,632.79	45	13	1,542,911.92	46	10
8	1,173,466.25	43	26	1,623,962.93	46	10
9	1,015,405.99	46	10	1,485,155.01	46	10
10	1,137,581.98	43	25	1,591,569.03	46	10

. In these solutions, the number of single-flight routes represents the number of dead flights, which incur additional costs for the airline since the aircraft has to fly to the starting city of the other route without passengers. The solutions obtained by considering the CA for the same problem are given in

Table 9. Results of proposed methods for 180 flights with codeshare flights.

	SA + MC			CS + MC
Run	Profit (USD)	Route Number	Codeshare Flights	Profit (USD)	Route Number	Codeshare Flights
1	1,641,346.89	38	52	1,616,784.05	44	24
2	1,561,180.42	36	62	1,885,131.11	42	28
3	1,605,142.62	36	62	1,807,812.28	43	30
4	1,626,210.07	36	54	1,912,243.16	45	14
5	1,731,738.91	37	56	1,825,287.15	40	37
6	1,659,333.69	39	50	1,947,504.96	46	12
7	1,548,412.45	36	54	1,812,634.21	43	33
8	1,571,768.08	37	55	1,815,737.61	43	26
9	1,591,618.95	36	57	1,819,234.24	44	20
10	1,501,147.29	35	64	1,870,894.85	42	28

, including the number of routes and the number of flights made with the CA. The effect of codesharing on the profits in each algorithm is shown in Figure 5a,b. These results confirm that codesharing has a significant contribution to increasing the airline’s profit.

Since airline operations have a complex structure, there are differences between the planned and actual schedule. Evaluating scenarios containing uncertainty can minimize these differences. The results obtained in this study were obtained by using many scenarios and two algorithms. Looking at the literature, software programs are available to simulate airline operations. For example, open data and open source simulation programs are available for air traffic management (AirTrafficSim, BlueSky). The programs quickly simulate operations and take into account a large number of scenarios. They are easy-to-use and accessible programs [78,79]. The effectiveness of the results in the scenarios we produced in our study can be seen by testing them with such simulators. Additionally, better results can be achieved with more scenarios.

6. Conclusions

The aviation industry has witnessed a steady increase in demand owing to its speed, convenience, and reliability. Airlines face numerous challenges to serve a vast network of destinations. Foremost among these, there are the high costs of investing in a large fleet, personnel, and equipment. To mitigate these costs, airlines have resorted to forming strategic partnerships through codesharing, enabling them to leverage each other’s flights.

In addition to this, airlines face the challenge of optimizing their operations across multiple domains such as FS, FA, and AR. While solving each of these problems independently may seem like an option, doing so may lead to local optima and increase operational costs. To overcome this, it is imperative to simultaneously address all these problems in a model that yields a more effective solution.

To address these challenges, this study proposes a novel model that integrates the FS, FA, and AR problems with codesharing options. Specifically, a three-stage stochastic non-linear mixed-integer mathematical model is developed. The first step involves determining the flights that will be operated by the airline and those that will be operated under the CA. Based on this, an appropriate fleet type is assigned to each flight, and the routes are created. In terms of the second stage, the estimated departure time for the airline’s own flights is selected from a predetermined time window. In the final step, the timetable of the routes assigned to each aircraft is updated, and the number of passengers to be carried by the aircraft is determined.

However, another critical challenge that makes solving the aforementioned optimization problems complex is the uncertainty surrounding customer demand and non-cruise times. Capturing these uncertainties accurately is crucial in modeling and decision-making problems. The developed model has the ability to model uncertainties in demand and non-cruise time parameters stochastically. The variation in the values of stochastic parameters will have an impact on the decisions to be made. Therefore, evaluating a wide range of scenarios reflecting the different values of stochastic parameters will ensure robust results. In this study, the Monte Carlo technique is used to analyze a vast number of scenarios for stochastic parameters.

On the other hand, the increase in the number of flights and scenarios makes it impractical to solve the proposed model using optimization software. Therefore, in this study, the SA and CS algorithms are proposed for solving the non-linear mixed-integer model by integrating the Monte Carlo technique. Testing these algorithms on real-life data reveals that the proposed model accurately represents the problem and produces efficient and robust results with the proposed solution algorithms. The CS algorithm is shown to have a better performance than SA for the integrated FS, FA, and AR problems.

The proposed model assumes that all fleet types can use all airports. However, in reality, certain fleets may not be able to use certain airports for technical reasons, and this restriction is referred to as station purity. This limitation serves as a constraint for the proposed model. Additionally, the developed model neglects the option of overbooking, where airlines sell more tickets than the airplane’s capacity. The model considers a single fare class but does not account for fare classes such as business class and first class. All these limitations present research potential for future studies. Finally, techniques such as Sample Average Approximation could be implemented to the metaheuristics to analyze stochastic parameters.