The multi-commodity flow problem

JuMP implementation of the multicommodity transportation model AMPL: A Modeling Language for Mathematical Programming, 2nd ed by Robert Fourer, David Gay, and Brian W. Kernighan 4-1.

Originally contributed by Louis Luangkesorn, February 26, 2015.

using JuMP
import HiGHS
import Test

function example_multi(; verbose = true)
    orig = ["GARY", "CLEV", "PITT"]
    dest = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
    prod = ["bands", "coils", "plate"]
    numorig = length(orig)
    numdest = length(dest)
    numprod = length(prod)
    # supply(prod, orig) amounts available at origins
    supply = [
        400 700 800
        800 1600 1800
        200 300 300
    ]
    # demand(prod, dest) amounts required at destinations
    demand = [
        300 300 100 75 650 225 250
        500 750 400 250 950 850 500
        100 100 0 50 200 100 250
    ]
    # limit(orig, dest) of total units from any origin to destination
    limit = [625.0 for j in 1:numorig, i in 1:numdest]
    # cost(dest, orig, prod) Shipment cost per unit
    cost = reshape(
        [
            [
                [30, 10, 8, 10, 11, 71, 6]
                [22, 7, 10, 7, 21, 82, 13]
                [19, 11, 12, 10, 25, 83, 15]
            ]
            [
                [39, 14, 11, 14, 16, 82, 8]
                [27, 9, 12, 9, 26, 95, 17]
                [24, 14, 17, 13, 28, 99, 20]
            ]
            [
                [41, 15, 12, 16, 17, 86, 8]
                [29, 9, 13, 9, 28, 99, 18]
                [26, 14, 17, 13, 31, 104, 20]
            ]
        ],
        7,
        3,
        3,
    )
    # DECLARE MODEL
    multi = Model(HiGHS.Optimizer)
    # VARIABLES
    @variable(multi, trans[1:numorig, 1:numdest, 1:numprod] >= 0)
    # OBJECTIVE
    @objective(
        multi,
        Max,
        sum(
            cost[j, i, p] * trans[i, j, p] for i in 1:numorig, j in 1:numdest,
            p in 1:numprod
        )
    )
    # CONSTRAINTS
    # Supply constraint
    @constraint(
        multi,
        supply_con[i in 1:numorig, p in 1:numprod],
        sum(trans[i, j, p] for j in 1:numdest) == supply[p, i]
    )
    # Demand constraint
    @constraint(
        multi,
        demand_con[j in 1:numdest, p in 1:numprod],
        sum(trans[i, j, p] for i in 1:numorig) == demand[p, j]
    )
    # Total shipment constraint
    @constraint(
        multi,
        total_con[i in 1:numorig, j in 1:numdest],
        sum(trans[i, j, p] for p in 1:numprod) - limit[i, j] <= 0
    )
    optimize!(multi)
    Test.@test termination_status(multi) == OPTIMAL
    Test.@test primal_status(multi) == FEASIBLE_POINT
    Test.@test objective_value(multi) == 225_700.0
    if verbose
        println("RESULTS:")
        for i in 1:length(orig)
            for j in 1:length(dest)
                for p in 1:length(prod)
                    print(
                        " $(prod[p]) $(orig[i]) $(dest[j]) = $(value(trans[i, j, p]))\t",
                    )
                end
                println()
            end
        end
    end
    return
end

example_multi()

Presolving model
44 rows, 60 cols, 165 nonzeros
41 rows, 60 cols, 149 nonzeros
Presolve : Reductions: rows 41(-10); columns 60(-3); elements 149(-40)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.6599957764e+03 Ph1: 41(149); Du: 60(1660) 0s
         60    -2.2570000000e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 60
Objective value     :  2.2570000000e+05
HiGHS run time      :          0.00
RESULTS:
 bands GARY FRA = 25.0	 coils GARY FRA = 500.0	 plate GARY FRA = 100.0
 bands GARY DET = 125.0	 coils GARY DET = 0.0	 plate GARY DET = 50.0
 bands GARY LAN = 0.0	 coils GARY LAN = 0.0	 plate GARY LAN = 0.0
 bands GARY WIN = 0.0	 coils GARY WIN = 0.0	 plate GARY WIN = 50.0
 bands GARY STL = 250.0	 coils GARY STL = 300.0	 plate GARY STL = 0.0
 bands GARY FRE = 0.0	 coils GARY FRE = 0.0	 plate GARY FRE = 0.0
 bands GARY LAF = 0.0	 coils GARY LAF = 0.0	 plate GARY LAF = 0.0
 bands CLEV FRA = 275.0	 coils CLEV FRA = 0.0	 plate CLEV FRA = 0.0
 bands CLEV DET = 0.0	 coils CLEV DET = 300.0	 plate CLEV DET = 50.0
 bands CLEV LAN = 100.0	 coils CLEV LAN = 0.0	 plate CLEV LAN = -0.0
 bands CLEV WIN = 0.0	 coils CLEV WIN = 0.0	 plate CLEV WIN = 0.0
 bands CLEV STL = 0.0	 coils CLEV STL = 625.0	 plate CLEV STL = 0.0
 bands CLEV FRE = 225.0	 coils CLEV FRE = 400.0	 plate CLEV FRE = 0.0
 bands CLEV LAF = 100.0	 coils CLEV LAF = 275.0	 plate CLEV LAF = 250.0
 bands PITT FRA = 0.0	 coils PITT FRA = 0.0	 plate PITT FRA = 0.0
 bands PITT DET = 175.0	 coils PITT DET = 450.0	 plate PITT DET = -0.0
 bands PITT LAN = 0.0	 coils PITT LAN = 400.0	 plate PITT LAN = 0.0
 bands PITT WIN = 75.0	 coils PITT WIN = 250.0	 plate PITT WIN = 0.0
 bands PITT STL = 400.0	 coils PITT STL = 25.0	 plate PITT STL = 200.0
 bands PITT FRE = 0.0	 coils PITT FRE = 450.0	 plate PITT FRE = 100.0
 bands PITT LAF = 150.0	 coils PITT LAF = 225.0	 plate PITT LAF = 0.0