multi_sets.jl

#############################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# Author: Louis Luangkesorn
# Date: February 26, 2015
#
# 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 multi.mod and multi.dat
#############################################################################

using Compat, JuMP, Clp, Cbc, GLPKMathProgInterface, DataArrays, DataFrames


function MultiMod()

    function PrintSolution(status, Trans, ORIG, DEST, PROD)
        println("\nRESULTS:")
        separator = Dict(:bands => "\t", :coils => "\t", :plate => "\n")
        if status == :Optimal
            println("\nObjective: $(getobjectivevalue(multi))")
            println("\nVariables:")
            for o in ORIG, d in DEST, p in PROD
                print("$p $o $d = $(getvalue(Trans[p,o,d]))$(separator[p])")
            end
            println("\nShadow prices of demand:")
            for d in DEST, p in PROD
                @printf("%s %s = %5.1f%s", p, d, getdual(c_demand[p,d]), separator[p])
            end
        else
            println("  No solution")
        end
        println("")
    end

    function PrintSolutionDataFrame(status, Trans, ORIG, DEST, PROD)
        trans = DataFrame(prod = Symbol[], orig = Symbol[], dest = Symbol[], value = Float64[])
        for p in PROD, o in ORIG, d in DEST
            push!(trans, @data([p,o,d,getvalue(Trans[p,o,d])]))
        end
        trans = unstack(trans, :dest, :value)

        show(IOContext(STDOUT, displaysize=(100,100)), "text/plain", trans)
    end

    # index sets
    orig = [:GARY, :CLEV, :PITT]
    dest = [:FRA, :DET, :LAN, :WIN, :STL, :FRE, :LAF]
    prod = [:bands, :coils, :plate]

    # Read input data into DataFrames and then into Dictionaries with 2D or 3D tuples of index sets as keys
    # wsv"""...""" makes a DataFrame of the multiline string

    # supply(prod, orig) amounts available at origins
    supplytable = wsv"""
    prod    GARY    CLEV    PITT
    bands   400     700     800
    coils   800     1600    1800
    plate   200     300     300
    """
    supply = Dict((Symbol(r[:prod]),o) => r[o] for r in eachrow(supplytable), o in orig)

    # demand(prod, dest) amounts required at destinations
    demandtable = wsv"""
    prod    FRA DET LAN WIN STL FRE LAF
    bands   300 300 100 75  650 225 250
    coils   500 750 400 250 950 850 500
    plate   100 100 0   50  200 100 250
    """
    demand = Dict((Symbol(r[:prod]),d) => r[d] for r in eachrow(demandtable), d in dest)

    # limit(orig, dest) of total units from any origin to destination
    defaultlimit = float(625)
    limit = Dict((o,d) => defaultlimit::AbstractFloat for o in orig, d in dest)

    # cost(prod, orig, dest) Shipment cost per unit
    costtable = wsv"""
    prod    orig    FRA DET LAN WIN STL FRE LAF
    bands   GARY    30  10  8   10  11  71  6
    bands   CLEV    22  7   10  7   21  82  13
    bands   PITT    19  11  12  10  25  83  15

    coils   GARY    39  14  11  14  16  82  8
    coils   CLEV    27  9   12  9   26  95  17
    coils   PITT    24  14  17  13  28  99  20

    plate   GARY    41  15  12  16  17  86  8
    plate   CLEV    29  9   13  9   28  99  18
    plate   PITT    26  14  17  13  31  104 20
    """
    cost = Dict((Symbol(r[:prod]),Symbol(r[:orig]),d) => r[d] for r in eachrow(costtable), d in dest)

    multi = Model(solver=ClpSolver())   # or CbcSolver() or GLPKSolverLP()

    @variables multi begin
        Trans[p in prod, o in orig, d in dest] >= 0
    end

    @constraints multi begin
        c_supply[p in prod, o in orig],
            sum(Trans[p,o,d] for d in dest) == supply[p,o]

        c_demand[p in prod, d in dest],
            sum(Trans[p,o,d] for o in orig) == demand[p,d]

        c_total_shipment[o in orig, d in dest],
            sum(Trans[p,o,d] for p in prod) - limit[o,d] <= 0
    end

    @objective multi Max begin
        sum(cost[p,o,d] * Trans[p,o,d] for p in prod, o in orig, d in dest)
    end

    status = solve(multi)
    if status == :Optimal
        result = Trans
        PrintSolution(status, Trans, orig, dest, prod)
        #PrintSolutionDataFrame(status, Trans, orig, dest, prod)
    else
        result = status
    end

end

MultiMod()