Airline seat allocation problem#
Attribution#
The following problem statement is adapted from an exercise and examples presented by Birge and Louveaux (2011).
Birge, J. R., & Louveaux, F. (2011). Introduction to stochastic programming. Springer Science & Business Media.
The adaptations include a change to parameters for consistency among reformulations of the problem, and additional treatments for sample average approximation (SAA) with chance constraints.
Problem description#
An airline is deciding how to partition a new plane for the AmsterdamBuenos Aires route. This plane can seat 200 economyclass passengers. A section can be created for firstclass seats, but each of these seats takes the space of 2 economyclass seats. A businessclass section can also be created, but each of these takes the space of 1.5 economyclass seats. The profit for a firstclass ticket is three times the profit of an economy ticket, while a businessclass ticket has a profit of two times an economy ticketās profit. Once the plane is partitioned into these seating classes, it cannot be changed.
The airlines knows that the plane will not always be full in every section. The airline has initially identified three scenarios to consider with about equal frequency:
Weekday morning and evening traffic;
Weekend traffic; and
Weekday midday traffic.
Under Scenario 1 the airline thinks they can sell as many as 20 firstclass tickets, 50 businessclass tickets, and 200 economy tickets. Under Scenario 2 these figures are 10 , 24, and 175, while under Scenario 3, they are 6, 10, and 150, respectively. The following table summarizes the forecast demand for these three scenarios.
Scenario 
Firstclass seats 
Businessclass seats 
Economyclass seats 

(1) weekday morning and evening 
20 
50 
200 
(2) weekend 
10 
24 
175 
(3) weekday midday 
6 
10 
150 
Average Scenario 
12 
28 
175 
The goal of the airline is to maximize ticket revenue. For marketing purposes, the airline will not sell more tickets than seats in each of the sections (hence no overbooking strategy). We further assume customers seeking a firstclass or businessclass seat will not downgrade if those seats are unavailable.
Installation and imports#
# install dependencies and select solver
%pip install q amplpy pandas matplotlib numpy scipy
SOLVER = "highs"
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["highs"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
Using default Community Edition License for Colab. Get yours at: https://ampl.com/ce
Licensed to AMPL Community Edition License for the AMPL Model Colaboratory (https://colab.ampl.com).
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Problem data#
Pandas DataFrames and Series are used to encode problem data in the following cell, and to encode problem solutions in subsequent cells.
# scenario data
demand = pd.DataFrame(
{
"morning and evening": {"F": 20, "B": 50, "E": 200},
"weekend": {"F": 10, "B": 24, "E": 175},
"midday": {"F": 6, "B": 10, "E": 150},
}
).T
# global revenue and seat factor data
capacity = 200
revenue_factor = pd.Series({"F": 3.0, "B": 2.0, "E": 1.0})
seat_factor = pd.Series({"F": 2.0, "B": 1.5, "E": 1.0})
Analytics#
Prior to optimization, a useful first step is to prepare an analytics function to display performance for any given allocation of seats. The firststage decision variables are the number of seats allocated for each class \(c\in C\). We would like to provide a analysis showing the operational consequences for any proposed allocation of seats. For this purpose, we create a function seat_report()
that show the tickets that can be sold in each scenario, the resulting revenue, and the unsatisfied demand (āspillageā).
To establish a basis for analyzing possible solutions to the airlineās problem, this function first is demonstrated for the case where the airplane is configured as entirely economyclass.
# function to report analytics for any given seat allocations
def seat_report(seats, demand):
classes = seats.index
# report seat allocation
equivalent_seats = pd.DataFrame(
{
"seat allocation": {c: seats[c] for c in classes},
"economy equivalent seat allocation": {
c: seats[c] * seat_factor[c] for c in classes
},
}
).T
equivalent_seats["TOTAL"] = equivalent_seats.sum(axis=1)
print("\nSeat Allocation")
display(equivalent_seats)
# tickets sold is the minimum of available seats and demand
tickets = pd.DataFrame()
for c in classes:
tickets[c] = np.minimum(seats[c], demand[c])
print("\nTickets Sold")
display(tickets)
# seats unsold
unsold = pd.DataFrame()
for c in classes:
unsold[c] = seats[c]  tickets[c]
print("\nSeats not Sold")
display(unsold)
# spillage (unmet demand)
spillage = demand  tickets
print("\nSpillage (Unfulfilled Demand)")
display(spillage)
# compute revenue
revenue = tickets.dot(revenue_factor)
print(
f"\nExpected Revenue (in units of economy ticket price): {revenue.mean():.2f}"
)
# charts
fig, ax = plt.subplots(2, 1, figsize=(8, 6))
revenue.plot(ax=ax[0], kind="barh", title="Revenue by Scenario")
ax[0].plot([revenue.mean()] * 2, ax[0].get_ylim(), "", lw=1.4)
ax[0].set_xlabel("Revenue")
tickets[classes].plot(
ax=ax[1], kind="barh", rot=0, stacked=False, title="Demand Scenarios"
)
demand[classes].plot(
ax=ax[1],
kind="barh",
rot=0,
stacked=False,
title="Demand Scenarios",
alpha=0.2,
)
for c in classes:
ax[1].plot([seats[c]] * 2, ax[1].get_ylim(), "", lw=1.4)
ax[1].set_xlabel("Seats")
fig.tight_layout()
return
# a trial solution
seats_all_economy = pd.Series({"F": 0, "B": 0, "E": 200})
seat_report(seats_all_economy, demand)
Seat Allocation
F  B  E  TOTAL  

seat allocation  0.0  0.0  200.0  200.0 
economy equivalent seat allocation  0.0  0.0  200.0  200.0 
Tickets Sold
F  B  E  

morning and evening  0  0  200 
weekend  0  0  175 
midday  0  0  150 
Seats not Sold
F  B  E  

morning and evening  0  0  0 
weekend  0  0  25 
midday  0  0  50 
Spillage (Unfulfilled Demand)
F  B  E  

morning and evening  20  50  0 
weekend  10  24  0 
midday  6  10  0 
Expected Revenue (in units of economy ticket price): 175.00
Model 1. Deterministic solution for the average demand scenario#
A common starting point in stochastic optimization is to solve the deterministic problem where future demands are fixed at their mean values and compute the corresponding optimal solution. The resulting value of the objective has been called the expectation of the expected value problem (EEV) by Birge, or the expected value of the mean (EVM) solution by others.
Let us introduce the set \(C\) of the three possible classes, i.e., \(C=\{F,B,E\}\). The objective function is to maximize ticket revenue.
where \(r_c\) is the revenue from selling a ticket for a seat in class \(c\in C\).
Let \(s_c\) denote the number of seats of class \(c \in C\) installed in the new plane. Let \(f_c\) be the scale factor denoting the number of economy seats displaced by one seat in class \(c \in C\). Then, since there is a total of 200 economyclass seats that could fit on the plane, the capacity constraint reads as
Let \(\mu_c\) be the mean demand for seats of class \(c\in C\), and let \(t_c\) be the number of tickets of class \(c\in C\) that are sold. To ensure we do not sell more tickets than available seats nor more than demand, we need to add two more constraints:
Lastly, both ticket and seat variables need to be nonnegative integers, so we add the constraints \(\bm{t}, \bm{s} \in \mathbb{Z}_+\).
The following cell presents an AMPL model implementing this model.
%%writefile airline_deterministic.mod
param capacity;
set CLASSES;
param seat_factor{CLASSES};
param demand{CLASSES};
param revenue_factor{CLASSES};
# first stage variables and constraints
var seats{CLASSES} integer >= 0;
s.t. plane_seats: sum{c in CLASSES}(seats[c] * seat_factor[c]) <= capacity;
# second stage variable and constraints
var tickets{CLASSES} integer >= 0;
s.t. demand_limits {c in CLASSES}: tickets[c] <= demand[c];
s.t. seat_limits {c in CLASSES}: tickets[c] <= seats[c];
# objective
maximize revenue: sum{c in CLASSES}(tickets[c] * revenue_factor[c]);
Overwriting airline_deterministic.mod
def airline_deterministic(demand):
# Create AMPL instance and load the model
ampl = AMPL()
ampl.read("airline_deterministic.mod")
# load the data
ampl.set["CLASSES"] = demand.columns.tolist()
ampl.param["demand"] = demand.mean()
ampl.param["revenue_factor"] = revenue_factor
ampl.param["seat_factor"] = seat_factor
ampl.param["capacity"] = capacity
# set solver
ampl.option["solver"] = SOLVER
return ampl
# Solve a given model and return a Pandas series of the seats for each class
def airline_solve(model):
model.solve()
return pd.Series(model.var["seats"].to_dict()).reindex(index=["F", "B", "E"])
# Solve deterministic model to obtain the expectation of the expected value problem (EEV)
model_eev = airline_deterministic(demand)
seats_eev = airline_solve(model_eev)
seat_report(seats_eev, demand)
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 226
2 simplex iterations
1 branching nodes
Seat Allocation
F  B  E  TOTAL  

seat allocation  12.0  28.0  134.0  174.0 
economy equivalent seat allocation  24.0  42.0  134.0  200.0 
Tickets Sold
F  B  E  

morning and evening  12  28  134 
weekend  10  24  134 
midday  6  10  134 
Seats not Sold
F  B  E  

morning and evening  0  0  0 
weekend  2  4  0 
midday  6  18  0 
Spillage (Unfulfilled Demand)
F  B  E  

morning and evening  8  22  66 
weekend  0  0  41 
midday  0  0  16 
Expected Revenue (in units of economy ticket price): 203.33
Model 2. Twostage stochastic optimization and its extensive form#
If we assume demand is not certain, we can formulate a twostage stochastic optimization problem. The firststage or hereandnow variables are the \(s_c\)ās, those related to the seat allocations. Due to their dependence on the realized demand \(\boldsymbol{z}\), the variables \(t_c\)ās describing the number of tickets sold, are secondstage or recourse decision variables. The full problem formulation is as follows: the first stage problem is
where \(Q(\boldsymbol{s},\boldsymbol{z})\) is the value of the secondstage problem, defined as
In view of the assumption that there is only a finite number \(N=S\) of scenarios for ticket demand, we can write the extensive form of the twostage stochastic optimization problem above and solve it exactly. To do so, we modify the secondstage variables \(t_{c,s}\) so that they are indexed by both class \(c\) and scenario \(s\). The expectation can thus be replaced with the average revenue over the \(N\) scenarios, that is
where the fraction \(\frac{1}{N}\) appears since we assume that all \(N\) scenarios are equally likely. The first stage constraint remains unchanged, while the second stage constraints are tuplicated for each scenario \(s \in S\), namely
The following cell presents an AMPL model implementing this model and solving it for the \(N=3\) equiprobable scenarios introduced above.
%%writefile airline_stochastic.mod
param capacity;
set CLASSES;
set SCENARIOS;
param demand{CLASSES, SCENARIOS};
param seat_factor{CLASSES};
param revenue_factor{CLASSES};
# first stage variables and constraints
var seats{CLASSES} integer >= 0;
s.t. plane_seats: sum{c in CLASSES}(seats[c] * seat_factor[c]) <= capacity;
# second stage variable and constraints
var tickets{CLASSES, SCENARIOS} integer >= 0;
s.t. demand_limits {c in CLASSES, s in SCENARIOS}: tickets[c, s] <= demand[c, s];
s.t. seat_limits {c in CLASSES, s in SCENARIOS}: tickets[c, s] <= seats[c];
# objective
maximize revenue: sum{c in CLASSES, s in SCENARIOS}(tickets[c, s] * revenue_factor[c]);
Overwriting airline_stochastic.mod
def airline_stochastic(demand):
# Create AMPL instance and load the model
ampl = AMPL()
ampl.read("airline_stochastic.mod")
# load the data
ampl.set["CLASSES"] = demand.columns.tolist()
ampl.set["SCENARIOS"] = demand.index.values.tolist()
ampl.param["demand"] = demand.T
ampl.param["revenue_factor"] = revenue_factor
ampl.param["seat_factor"] = seat_factor
ampl.param["capacity"] = capacity
# set solver
ampl.option["solver"] = SOLVER
return ampl
# create and solve model for the three scenarios defined above
model_stochastic = airline_stochastic(demand)
seats_stochastic = airline_solve(model_stochastic)
seat_report(seats_stochastic, demand)
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 628
7 simplex iterations
1 branching nodes
Seat Allocation
F  B  E  TOTAL  

seat allocation  10.0  20.0  150.0  180.0 
economy equivalent seat allocation  20.0  30.0  150.0  200.0 
Tickets Sold
F  B  E  

morning and evening  10  20  150 
weekend  10  20  150 
midday  6  10  150 
Seats not Sold
F  B  E  

morning and evening  0  0  0 
weekend  0  0  0 
midday  4  10  0 
Spillage (Unfulfilled Demand)
F  B  E  

morning and evening  10  30  50 
weekend  0  4  25 
midday  0  0  0 
Expected Revenue (in units of economy ticket price): 209.33
Model 3. Adding chance constraints#
The airline wishes a special guarantee for its clients enrolled in its loyalty program. In particular, it wants \(98\%\) probability to cover the demand of firstclass seats and \(95\%\) probability to cover the demand of businessclass seats (by clients of the loyalty program). Firstclass passengers are covered if they can purchase a firstclass seat. Businessclass passengers are covered if they purchase a businessclass seat or upgrade to a firstclass seat.
Assume the demand of loyaltyprogram passengers is normally distributed as \(z_F \sim \mathcal N(\mu_F, \sigma_F^2)\) and \(z_B \sim \mathcal N(\mu_B, \sigma_B^2)\) for firstclass and business, respectively, where the parameters are given in the table below. For completeness, we also include the parameters for economyclass passengers.
\(\mu_{\cdot}\) 
\(\sigma_{\cdot}\) 


F 
12 
4 
B 
28 
8 
E 
175 
20 
We further assume that the demands for firstclass and businessclass seats are independent of each other and of the scenario (time of the week).
Let \(s_F\) be the number of firstclass seats and \(s_B\) the number of business seats. The probabilistic constraints are
These are can be written equivalently as linear constraints, specifically
For the second constraint, we use the fact that the sum of the two independent random variables \(z_F\) and \(z_B\) is normally distributed with mean \(\mu_F + \mu_B\) and variance \(\sigma_F^2 + \sigma_B^2\).
We add these equivalent linear counterparts of the two chance constraints to the stochastic optimization model. Rather than writing a function to create a whole new model, we can use the prior function to create and add the two chance constraints using decorators.
# import the scipy.stats.norm object needed to use the quantile function of the standard normal distribution
from scipy.stats import norm
# define the mean and standard deviation of the demand for each class
mu = demand.mean()
sigma = {"F": 4, "B": 16, "E": 20}
display(pd.DataFrame({"mu": mu, "sigma": sigma}))
# create a new model with chance constraints that takes as input also
# the target quality of service (QoS) levels for classes F and B
def airline_cc(demand, QoSF=0.98, QoSFB=0.95):
# create twostage stochastic model as before
m = airline_stochastic(demand)
# add equivalent counterparts of the two chance constraints to the first stage problem
# the two coefficients related the inverse CDF of the standard normal are computed using the norm.ppf function
first_class = "s.t. first_class: seats['F']  {} >= {};".format(
mu["F"], norm.ppf(QoSF) * sigma["F"]
)
m.eval(first_class)
business_class = "s.t. business_class: seats['F'] + seats['B']  {} >= {};".format(
mu["F"] + mu["B"], norm.ppf(QoSFB) * np.sqrt(sigma["F"] ** 2 + sigma["B"] ** 2)
)
m.eval(business_class)
print(m.con["first_class"])
print(m.con["business_class"])
return m
# create and solve model
model_cc = airline_cc(demand)
seats_cc = airline_solve(model_cc)
seat_report(seats_cc, demand)
mu  sigma  

F  12.0  4 
B  28.0  16 
E  175.0  20 
subject to first_class: seats['F']  12 >= 8.21499564252729;
subject to business_class: seats['F'] + seats['B']  40 >=
27.127620970404912;
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 531
4 simplex iterations
1 branching nodes
Seat Allocation
F  B  E  TOTAL  

seat allocation  21.0  47.0  87.0  155.0 
economy equivalent seat allocation  42.0  70.5  87.0  199.5 
Tickets Sold
F  B  E  

morning and evening  20  47  87 
weekend  10  24  87 
midday  6  10  87 
Seats not Sold
F  B  E  

morning and evening  1  0  0 
weekend  11  23  0 
midday  15  37  0 
Spillage (Unfulfilled Demand)
F  B  E  

morning and evening  0  3  113 
weekend  0  0  88 
midday  0  0  63 
Expected Revenue (in units of economy ticket price): 177.00
Model 4. Solving the case of continuous demand distributions using the SAA method#
Let us now move past the simplifying assumption that there are only three equally likely scenarios and consider the case where the demand is described by a random vector \((z_F, z_B, z_E)\), where \(z_c\) is the demand for seats of class \(c\in C\). The demand for class \(c\) is assumed to be independent and normally distributed with mean \(\mu_c\) and variance \(\sigma_c^2\) as reported in the following table
\(\mu\) 
\(\sigma\) 


F 
12 
4 
B 
28 
8 
E 
175 
20 
Note that we model the demand for each class using a continuous random variable, which is a simplification of the real world, where the ticket demand is always a discrete nonnegative number. Therefore, we round down all the obtained random numbers.
However, now that the number of scenarios is not finite anymore, we cannot solve the problem exactly. Instead, we can use the SAA method to approximate the expected value appearing in the objective function. The first step of the SAA is to generate a collection of \(N\) scenarios \((z_{F,s}, z_{B,s}, z_{E,s})\) for \(s=1,\ldots,N\). We can do this by sampling from the normal distributions with the given means and variances.
For sake of generality, we create a script to generate scenarios in which the three demands have a general correlation structure captured by a correlation matrix \(\bm{\rho}\).
Model 5. Adding correlations between different demand types#
Now assume the ticket demand for the three categories is captured by a \(3\)dimensional multivariate normal distribution mean \(\mu=(\mu_F, \mu_B, \mu_E)\), variances \((\sigma_F^2, \sigma_B^2, \sigma_E^2)\) and a symmetric correlation matrix
The covariance matrix is given by \(\Sigma = \text{diag}(\sigma)\ P\ \text{diag}(\sigma)\) or
We assume \(\rho_{FB} = 0.6\), \(\rho_{BE} = 0.4\) and \(\rho_{FE} = 0.2\).
We now sample \(N=1000\) scenarios from this multivariate correlated normal distribution and use the SAA method to approximate the solution of the twostage stochastic optimization problem.
# sample size
N = 1000
# define the mean mu and standard deviation sigma of the demand for each class
mu = demand.mean()
sigma = {"F": 4, "B": 16, "E": 20}
display(pd.DataFrame({"mu": mu, "sigma": sigma}))
# correlation matrix
P = np.array([[1, 0.6, 0.2], [0.6, 1, 0.4], [0.2, 0.4, 1]])
# build covariance matrix from covariances and correlations
s = np.array(list(sigma.values()))
S = np.diag(s) @ P @ np.diag(s)
# create samples
np.random.seed(1)
samples = np.random.multivariate_normal(list(mu), S, N).round()
# truncate to integers and nonnegative values
classes = demand.columns
demand_saa = pd.DataFrame(samples, columns=classes)
demand_saa[demand_saa < 0] = 0
df = pd.DataFrame(mu, columns=["mu"])
df["sample means"] = demand_saa.mean()
display(df)
df = pd.DataFrame(pd.Series(sigma), columns=["sigma"])
df["sample std dev"] = demand_saa.std()
display(df)
print("\nModel Covariance")
df = pd.DataFrame(S, index=classes, columns=classes)
display(df)
print("\nSample Covariance")
display(pd.DataFrame(demand_saa.cov()))
fig, ax = plt.subplots(3, 3, figsize=(10, 10))
for i, ci in enumerate(classes):
for j, cj in enumerate(classes):
if i == j:
demand_saa[ci].hist(ax=ax[i, i], bins=20)
ax[i, i].set_title(f"Histogram {ci}")
else:
ax[i, j].plot(demand_saa[ci], demand_saa[cj], ".", alpha=0.4)
ax[i, j].set_xlabel(ci)
ax[i, j].set_ylabel(cj)
ax[i, j].set_title(f"Convariance: {ci} vs {cj}")
fig.tight_layout()
mu  sigma  

F  12.0  4 
B  28.0  16 
E  175.0  20 
mu  sample means  

F  12.0  12.058 
B  28.0  28.355 
E  175.0  175.436 
sigma  sample std dev  

F  4  4.046852 
B  16  15.879708 
E  20  19.503597 
Model Covariance
F  B  E  

F  16.0  38.4  16.0 
B  38.4  256.0  128.0 
E  16.0  128.0  400.0 
Sample Covariance
F  B  E  

F  16.377013  40.352763  15.402114 
B  40.352763  252.165140  121.395616 
E  15.402114  121.395616  380.390294 
model_ssa = airline_stochastic(demand_saa)
seats_saa = airline_solve(model_ssa)
seat_report_saa(seats_saa, demand_saa)
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 210891
2103 simplex iterations
1 branching nodes
Seat Allocation
F  B  E  TOTAL  

seat allocation  11.0  20.0  148.0  179.0 
economy equivalent seat allocation  22.0  30.0  148.0  200.0 
Mean Tickets Sold
F 9.868
B 17.021
E 147.245
dtype: float64
Mean Seats not Sold
F 1.132
B 2.979
E 0.755
dtype: float64
Mean Spillage (Unfulfilled Demand)
F 2.190
B 11.334
E 28.191
dtype: float64
Expected Revenue (in units of economy ticket price): 210.89