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
import TestFormulation
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
├ solver: HiGHS
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: noneOur 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] ≤ 600and 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] ≥ 300Finally, 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 : 1.68000e+03
Relative gap : 0.00000e+00
Dual objective value : 1.68000e+03
* Work counters
Solve time (sec) : 2.05994e-04
Simplex iterations : 3
Barrier iterations : 0
Node count : -1
What's the optimal shipment?
Test.@test is_solved_and_feasible(model)
for p in P, m in M
println(p, " => ", m, ": ", value(x[p, m]))
endSeattle => 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