The cannery problem

Original author: Louis Luangkesorn, January 30, 2015.

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

using JuMP
import HiGHS
import JSON

Formulation

The cannery problem assumes we are optimizing the shipment 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 distance from plant to market is $d_{p,m}$.

With a little effort, 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.0
 x[San-Diego,Chicago] + x[San-Diego,Topeka] + x[San-Diego,New-York] ≤ 600.0

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.0
 x[Seattle,Topeka] + x[San-Diego,Topeka] ≥ 300.0
 x[Seattle,New-York] + x[San-Diego,New-York] ≥ 300.0

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))

\[ 1.7 x_{Seattle,Chicago} + 1.8 x_{Seattle,Topeka} + 2.5 x_{Seattle,New-York} + 1.8 x_{San-Diego,Chicago} + 1.4 x_{San-Diego,Topeka} + 2.5 x_{San-Diego,New-York} \]

Solution

optimize!(model)

solution_summary(model)

* Solver : HiGHS

* Status
  Termination status : OPTIMAL
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Message from the solver:
  "kHighsModelStatusOptimal"

* Candidate solution
  Objective value      : 1.68000e+03
  Objective bound      : 0.00000e+00
  Dual objective value : 1.68000e+03

* Work counters
  Solve time (sec)   : 3.66211e-04
  Simplex iterations : 3
  Barrier iterations : 0

What's the optimal shipment?

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

RESULTS:
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