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

License: Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

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 [ ]: