JuliaOpt

Optimization packages for the Julia language.


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 [ ]:
# First we need to load the JuMP package
using JuMP

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 [ ]:
nrm = 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 [ ]:
@defVar(nrm, 0 <= BOStoMDW_R <= 80)
@defVar(nrm, 0 <= BOStoMDW_D <= 120)
@defVar(nrm, 0 <= BOStoSFO_R <= 75)
@defVar(nrm, 0 <= BOStoSFO_D <= 100)
@defVar(nrm, 0 <= MDWtoSFO_R <= 60)
@defVar(nrm, 0 <= MDWtoSFO_D <= 110)
nrm

Great, now we are getting somewhere!

Objective

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

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

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 @addConstraint, which takes two arguments: the model, and an expression with an inequality in it: <=, ==, >=.

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

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

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 [ ]:
status = solve(nrm)
status  # Should be `:Optimal`

Inspect the solution

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

In [ ]:
getValue(BOStoMDW_R)
In [ ]:
getValue(BOStoMDW_D)
In [ ]:
getObjectiveValue(nrm)

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

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

In [ ]:
nrm2 = Model()
# 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.
@defVar(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!

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 [ ]:
# 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 (after the semicolon ;) to exclude all cases where
# the origin is equal to the destination
@setObjective(nrm2, Max, sum{prices[(o,d,c)]*x[o,d,c], 
    o in airports, d in airports, c in classes; o != d})
In [ ]:
# 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")
    @addConstraint(nrm2, 
        sum{x[o,d,c], o in airports, c in classes; o!=d} <= plane_cap)
end
nrm2
In [ ]:
# Now flights into the hub
for o in airports
    if o == :MDW
        continue
    end
    println("Adding constraint for $o to hub (MDW)")
    @addConstraint(nrm2, 
        sum{x[o,d,c], d in airports, c in classes; o!=d} <= plane_cap)
end
nrm2
In [ ]:
# Now just need to solve it
solve(nrm2)
In [ ]:
getValue(x)