The cannery problem

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial was originally contributed by Louis Luangkesorn.

This tutorial solves the cannery problem from Dantzig, Linear Programming and Extensions, Princeton University Press, Princeton, NJ, 1963. This class of problem is known as a transshipment problem.

The purpose of this tutorial is to demonstrate how to use JSON data in the formulation of a JuMP model.

Required packages

This tutorial requires the following packages:

using JuMP
import HiGHS
import JSON

Formulation

The cannery problem assumes we are optimizing the shipment of cases of cans from production plants $p \in P$ to markets $m \in M$.

Each production plant $p$ has a capacity $c_p$, and each market $m$ has a demand $d_m$. The shipping cost per case of cans from plant $p$ to market $m$ is $d_{p,m}$.

We wish to find the distribution plan $x_{p,m}$, the number of cases of cans to ship from plant $p$ to market $m$, for $p \in P$ and $m \in M$ that minimizes the shipping costs. We can formulate our problem as the following linear program:

\[\begin{aligned} \min & \sum\limits_{p \in P}\sum\limits_{m \in M} d_{p,m} x_{p,m} \\ \text{s.t.} & \sum\limits_{m \in M} x_{p,m} \le c_p, && \forall p \in P \\ & \sum\limits_{p \in P} x_{p,m} \ge d_m, && \forall m \in M \\ & x_{p,m} \ge 0, && \forall p \in P, m \in M \end{aligned}\]

Data

A key feature of the tutorial is to demonstrate how to load data from JSON.

For simplicity, we've hard-coded it below. But if the data was available as a .json file, we could use data = JSON.parsefile(filename) to read in the data.

data = JSON.parse("""
{
    "plants": {
        "Seattle": {"capacity": 350},
        "San-Diego": {"capacity": 600}
    },
    "markets": {
        "New-York": {"demand": 300},
        "Chicago": {"demand": 300},
        "Topeka": {"demand": 300}
    },
    "distances": {
        "Seattle => New-York": 2.5,
        "Seattle => Chicago": 1.7,
        "Seattle => Topeka": 1.8,
        "San-Diego => New-York": 2.5,
        "San-Diego => Chicago": 1.8,
        "San-Diego => Topeka": 1.4
    }
}
""")

Dict{String, Any} with 3 entries:
  "plants"    => Dict{String, Any}("Seattle"=>Dict{String, Any}("capacity"=>350…
  "distances" => Dict{String, Any}("San-Diego => New-York"=>2.5, "Seattle => To…
  "markets"   => Dict{String, Any}("Chicago"=>Dict{String, Any}("demand"=>300),…

Create the set of plants:

P = keys(data["plants"])

KeySet for a Dict{String, Any} with 2 entries. Keys:
  "Seattle"
  "San-Diego"

Create the set of markets:

M = keys(data["markets"])

KeySet for a Dict{String, Any} with 3 entries. Keys:
  "Chicago"
  "Topeka"
  "New-York"

We also need a function to compute the distance from plant to market:

distance(p::String, m::String) = data["distances"]["$(p) => $(m)"]

distance (generic function with 1 method)

JuMP formulation

Now we're ready to convert our mathematical formulation into a JuMP model.

First, create a new JuMP model. Since we have a linear program, we'll use HiGHS as our optimizer:

model = Model(HiGHS.Optimizer)

A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: HiGHS

Our decision variables are indexed over the set of plants and markets:

@variable(model, x[P, M] >= 0)

2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, ["Seattle", "San-Diego"]
    Dimension 2, ["Chicago", "Topeka", "New-York"]
And data, a 2×3 Matrix{VariableRef}:
 x[Seattle,Chicago]    x[Seattle,Topeka]    x[Seattle,New-York]
 x[San-Diego,Chicago]  x[San-Diego,Topeka]  x[San-Diego,New-York]

We need a constraint that each plant can ship no more than its capacity:

@constraint(model, [p in P], sum(x[p, :]) <= data["plants"][p]["capacity"])

1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, ["Seattle", "San-Diego"]
And data, a 2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
 x[Seattle,Chicago] + x[Seattle,Topeka] + x[Seattle,New-York] ≤ 350
 x[San-Diego,Chicago] + x[San-Diego,Topeka] + x[San-Diego,New-York] ≤ 600

and each market must receive at least its demand:

@constraint(model, [m in M], sum(x[:, m]) >= data["markets"][m]["demand"])

1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, ["Chicago", "Topeka", "New-York"]
And data, a 3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, ScalarShape}}:
 x[Seattle,Chicago] + x[San-Diego,Chicago] ≥ 300
 x[Seattle,Topeka] + x[San-Diego,Topeka] ≥ 300
 x[Seattle,New-York] + x[San-Diego,New-York] ≥ 300

Finally, our objective is to minimize the transportation distance:

@objective(model, Min, sum(distance(p, m) * x[p, m] for p in P, m in M));

Solution

Let's optimize and look at the solution:

optimize!(model)
solution_summary(model)

* Solver : HiGHS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "kHighsModelStatusOptimal"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 1.68000e+03
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : 1.68000e+03

* Work counters
  Solve time (sec)   : 3.99590e-04
  Simplex iterations : 3
  Barrier iterations : 0
  Node count         : -1

What's the optimal shipment?

for p in P, m in M
    println(p, " => ", m, ": ", value(x[p, m]))
end

Seattle => Chicago: 300.0
Seattle => Topeka: 0.0
Seattle => New-York: 0.0
San-Diego => Chicago: 0.0
San-Diego => Topeka: 300.0
San-Diego => New-York: 300.0