## Disclaimer¶

This notebook is only working under the versions:

• JuMP 0.19 (unreleased, but currently in master)

• MathOptInterface 0.4.1

• GLPK 0.6.0

Description: Shows how to solve a network revenue management problem using JuMP.

Author: Iain Dunning

# Airline Network Revenue Management¶

In the airline network revenue management problem we are trying to decide how many tickets for each origin-destination (O-D) pair to sell at each price level. The goal is to maximize revenue, and we cannot sell more tickets than there is demand for, or space on the planes for.

## Three Flight Problem¶

We'll start with a toy problem that has three origin-destination pairs, with two price classes for each pair. The three origin-destination pairs are BOS-MDW, MDW-SFO, or BOS-SFO via MDW. BOS stands for Boston, MDW is Chicago-Midway, and SFO is San Francisco. Each O-D pair has a "regular" and "discount" fare class. The data we will use is summarized as follows:

PLANE CAPACITY: 166

BOS-MDW
Regular  Discount
Price:  428      190
Demand: 80       120

BOS-SFO
Regular  Discount
Price:  642      224
Demand: 75       100

MDW-SFO
Regular  Discount
Price:  512      190
Demand: 60       110
In [1]:
push!(Base.LOAD_PATH,"D:\\MOI")

Out[1]:
3-element Array{Any,1}:
"C:\\Users\\joaquimgarcia\\AppData\\Local\\Julia-0.6.0\\local\\share\\julia\\site\\v0.6"
"C:\\Users\\joaquimgarcia\\AppData\\Local\\Julia-0.6.0\\share\\julia\\site\\v0.6"
"D:\\MOI"                                                                               
In [2]:
# Load JuMP
using JuMP
using MathOptInterface
# Load solver package
using GLPK
# shortcuts
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

Out[2]:
MathOptInterface.Utilities

### Model¶

Now we need to create a Model. A Model is an object that has associated variables and constraints. We can also pick and customize different solvers to use with the model. In this case we'll assume you a LP solver installed - JuMP will detect it and use it automatically.

In [3]:
nrm = Model(optimizer = GLPK.GLPKOptimizerLP())

Out[3]:
A JuMP Model

We can see that JuMP displays the model in a human-readable form.

### Variables¶

Now we can create our variables, in the optimization sense. Variables have a name, bounds, and type. For this problem, we need six continuous variables, each with a lower and upper bound on their value.

Here we will spell out each variable one-by-one - rest assured that later we will do something smarter that will scale to millions of variables!

In [4]:
@variables(nrm, begin
0 <= BOStoMDW_R <= 80
0 <= BOStoMDW_D <= 120
0 <= BOStoSFO_R <= 75
0 <= BOStoSFO_D <= 100
0 <= MDWtoSFO_R <= 60
0 <= MDWtoSFO_D <= 110
end)
nrm

Out[4]:
A JuMP Model

Great, now we are getting somewhere!

### Objective¶

The objective is to maximize the revenue. We set the objective with @objective, which takes three arguments: the model, the sense (Max or Min), and the expression.

In [5]:
@objective(nrm, Max, 428*BOStoMDW_R + 190*BOStoMDW_D +
642*BOStoSFO_R + 224*BOStoSFO_D +
512*MDWtoSFO_R + 190*MDWtoSFO_D)
nrm

Out[5]:
A JuMP Model

### Constraints¶

We have only two constraints, one per physical flight: that the number of people on each flight is less than the capacity of the planes.

We add constraints with @constraint, which takes two arguments: the model, and an expression with an inequality in it: <=, ==, >=.

(note that there are other forms of @constraint that can be useful sometimes - see the manual or other examples for details)

In [6]:
@constraint(nrm, BOStoMDW_R + BOStoMDW_D +
BOStoSFO_R + BOStoSFO_D <= 166)
@constraint(nrm, MDWtoSFO_R + MDWtoSFO_D +
BOStoSFO_R + BOStoSFO_D <= 166)
nrm

Out[6]:
A JuMP Model

Our model is complete!

## Solve¶

The easy part! Lets check out the finished model before solving it. We didn't specify a solver before, but JuMP knows we have an LP solver installed, so it will use that.

In [7]:
# solve problem
JuMP.optimize(nrm)

In [8]:
@show JuMP.hasresultvalues(nrm)
@show JuMP.terminationstatus(nrm) == MOI.Success
@show JuMP.primalstatus(nrm) == MOI.FeasiblePoint

JuMP.hasresultvalues(nrm) = true
JuMP.terminationstatus(nrm) == MOI.Success = true
JuMP.primalstatus(nrm) == MOI.FeasiblePoint = true

Out[8]:
true

## Inspect the solution¶

We can use getvalue() to inspect the value of solutions

In [9]:
JuMP.resultvalue(BOStoMDW_R)

Out[9]:
80.0
In [10]:
JuMP.resultvalue(BOStoMDW_D)

Out[10]:
11.0
In [11]:
JuMP.getobjectivevalue(nrm)

UndefVarError: getobjectivevalue not defined


## General Model¶

We'll now code our model in a more general way, to take any number of cities and flights. Hard-coding every variable would be painful and hard to update - it'd be better to index the variables, just like in mathematical notation.

Consider a generalized version of our first problem, where there are flights in both directions and one extra airport YYZ - Toronto!

Rather than hardcode data, we will generate some random data.

(If you don't understand all of this right away, thats OK - its not critical to understanding JuMP)

In [12]:
# Set the random seed to ensure we always
# get the same stream of 'random' numbers
srand(1988)

# Lets create a vector of symbols, one for each airport
airports = [:BOS, :MDW, :SFO, :YYZ]
num_airport = length(airports)

# We'll also create a vector of fare classes
classes = [:REG, :DIS]

# All the demand and price data for each triple of
# (origin, destination, class) will be stored in
# 'dictionaries', also known as 'maps'.
demand = Dict()
prices = Dict()

# Generate a demand and price for each pair of airports
# To keep the code simple we will generate info for
# nonsense flights like BOS-BOS and SFO-SFO, but they
# won't appear in our final model.
for origin in airports, dest in airports
# Generate demand:
#  - Regular demand is Uniform(50,90)
#  - Discount demand is Uniform(100,130)
demand[(origin,dest,:REG)] = rand(50:90)
demand[(origin,dest,:DIS)] = rand(100:130)
# Generate prices:
#  - Regular price is Uniform(400,700)
#  - Discount price is Uniform(150,300)
prices[(origin,dest,:REG)] = rand(400:700)
prices[(origin,dest,:DIS)] = rand(150:300)
end

# Finally set all places to have the same capacity
plane_cap = rand(150:200)

# Lets look at a sample demand at random
@show demand[(:BOS,:YYZ,:REG)]

demand[(:BOS, :YYZ, :REG)] = 90

Out[12]:
90

Now we have the data, we can build the model.

In [13]:
nrm2 = Model(optimizer = GLPK.GLPKOptimizerLP())
# Create a variable indexed by 3 things:
# * Origin
# * Destination
# * Class
# And set the upper bound for each variable differently.
# The origins and destinations should be indexed across
# the vector of airports, and the class should be indexed
# across the vector of classes.
# We'll draw the upper bounds from the demand dictionary.
@variable(nrm2, 0 <= x[o=airports,
d=airports,
c=classes] <= demand[(o,d,c)])
# Note that we don't have to split it up across multiple lines,
# its personal preference!

Out[13]:
3-dimensional JuMPArray{JuMP.VariableRef,3,...} with index sets:
Dimension 1, Symbol[:BOS, :MDW, :SFO, :YYZ]
Dimension 2, Symbol[:BOS, :MDW, :SFO, :YYZ]
Dimension 3, Symbol[:REG, :DIS]
And data, a 4×4×2 Array{JuMP.VariableRef,3}:
[:, :, :REG] =
x[BOS,BOS,REG]  x[BOS,MDW,REG]  x[BOS,SFO,REG]  x[BOS,YYZ,REG]
x[MDW,BOS,REG]  x[MDW,MDW,REG]  x[MDW,SFO,REG]  x[MDW,YYZ,REG]
x[SFO,BOS,REG]  x[SFO,MDW,REG]  x[SFO,SFO,REG]  x[SFO,YYZ,REG]
x[YYZ,BOS,REG]  x[YYZ,MDW,REG]  x[YYZ,SFO,REG]  x[YYZ,YYZ,REG]

[:, :, :DIS] =
x[BOS,BOS,DIS]  x[BOS,MDW,DIS]  x[BOS,SFO,DIS]  x[BOS,YYZ,DIS]
x[MDW,BOS,DIS]  x[MDW,MDW,DIS]  x[MDW,SFO,DIS]  x[MDW,YYZ,DIS]
x[SFO,BOS,DIS]  x[SFO,MDW,DIS]  x[SFO,SFO,DIS]  x[SFO,YYZ,DIS]
x[YYZ,BOS,DIS]  x[YYZ,MDW,DIS]  x[YYZ,SFO,DIS]  x[YYZ,YYZ,DIS]

JuMP displays the variable in a compact form - note that the upper bound is blank because the upper bound is different for each variable.

In [14]:
# The objective is just like before, except now we can use
# the sum() functionality provided by JuMP to sum over all
# combinations of origin/destination/class, and provide a
# filter to exclude all cases where
# the origin is equal to the destination
@objective(nrm2, Max, sum(prices[(o,d,c)]*x[o,d,c] for
o in airports, d in airports, c in classes if o != d))

In [15]:
# Next we'll add constraints on flights from the hub
# to any of the non-hub airports.
for d in airports
if d == :MDW
continue
end
println("Adding constraint for hub (MDW) to $d") @constraint(nrm2, sum(x[o,d,c] for o in airports, c in classes if o!=d) <= plane_cap) end nrm2  Adding constraint for hub (MDW) to BOS Adding constraint for hub (MDW) to SFO Adding constraint for hub (MDW) to YYZ  Out[15]: A JuMP Model In [16]: # Now flights into the hub for o in airports if o == :MDW continue end println("Adding constraint for$o to hub (MDW)")
@constraint(nrm2,
sum(x[o,d,c] for d in airports, c in classes if o!=d) <= plane_cap)
end
nrm2

Adding constraint for BOS to hub (MDW)
Adding constraint for SFO to hub (MDW)
Adding constraint for YYZ to hub (MDW)

Out[16]:
A JuMP Model
In [17]:
# solve problem
JuMP.optimize(nrm2)

In [18]:
JuMP.resultvalue.(x)

Out[18]:
3-dimensional JuMPArray{Float64,3,...} with index sets:
Dimension 1, Symbol[:BOS, :MDW, :SFO, :YYZ]
Dimension 2, Symbol[:BOS, :MDW, :SFO, :YYZ]
Dimension 3, Symbol[:REG, :DIS]
And data, a 4×4×2 Array{Float64,3}:
[:, :, :REG] =
0.0  55.0  46.0  81.0
86.0   0.0  75.0  63.0
54.0  90.0   0.0  38.0
0.0  62.0  61.0   0.0

[:, :, :DIS] =
0.0   0.0  0.0  0.0
42.0   0.0  0.0  0.0
0.0   0.0  0.0  0.0
0.0  59.0  0.0  0.0
In [ ]: