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:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
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.
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
push!(Base.LOAD_PATH,"D:\\MOI")
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"
# Load JuMP
using JuMP
using MathOptInterface
# Load solver package
using GLPK
# shortcuts
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
MathOptInterface.Utilities
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.
nrm = Model(optimizer = GLPK.GLPKOptimizerLP())
A JuMP Model
We can see that JuMP displays the model in a human-readable form.
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!
@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
A JuMP Model
Great, now we are getting somewhere!
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.
@objective(nrm, Max, 428*BOStoMDW_R + 190*BOStoMDW_D +
642*BOStoSFO_R + 224*BOStoSFO_D +
512*MDWtoSFO_R + 190*MDWtoSFO_D)
nrm
A JuMP Model
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)
@constraint(nrm, BOStoMDW_R + BOStoMDW_D +
BOStoSFO_R + BOStoSFO_D <= 166)
@constraint(nrm, MDWtoSFO_R + MDWtoSFO_D +
BOStoSFO_R + BOStoSFO_D <= 166)
nrm
A JuMP Model
Our model is complete!
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.
# solve problem
JuMP.optimize(nrm)
@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
true
We can use getvalue()
to inspect the value of solutions
JuMP.resultvalue(BOStoMDW_R)
80.0
JuMP.resultvalue(BOStoMDW_D)
11.0
JuMP.getobjectivevalue(nrm)
UndefVarError: getobjectivevalue not defined
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)
# 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
90
Now we have the data, we can build the model.
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!
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.
# 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))
# 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
A JuMP Model
# 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)
A JuMP Model
# solve problem
JuMP.optimize(nrm2)
JuMP.resultvalue.(x)
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