Mini MILP "benchmark"

using JuMP
import GLPK
import Cbc

function example_urban_plan(factory)
    model = Model(factory)
    set_silent(model)
    # x is indexed by row and column
    @variable(model, 0 <= x[1:5, 1:5] <= 1, Int)
    # y is indexed by R or C, the points, and an index in 1:5. Note how JuMP
    # allows indexing on arbitrary sets.
    rowcol = ["R", "C"]
    points = [5, 4, 3, -3, -4, -5]
    @variable(model, 0 <= y[rowcol, points, 1:5] <= 1, Int)
    # Objective - combine the positive and negative parts
    @objective(
        model,
        Max,
        sum(
            3 * (y["R", 3, i] + y["C", 3, i]) +
            1 * (y["R", 4, i] + y["C", 4, i]) +
            1 * (y["R", 5, i] + y["C", 5, i]) -
            3 * (y["R", -3, i] + y["C", -3, i]) -
            1 * (y["R", -4, i] + y["C", -4, i]) -
            1 * (y["R", -5, i] + y["C", -5, i]) for i in 1:5
        )
    )
    # Constrain the number of residential lots
    @constraint(model, sum(x) == 12)
    # Add the constraints that link the auxiliary y variables to the x variables
    for i in 1:5
        @constraints(model, begin
            # Rows
            y["R", 5, i] <= 1 / 5 * sum(x[i, :]) # sum = 5
            y["R", 4, i] <= 1 / 4 * sum(x[i, :]) # sum = 4
            y["R", 3, i] <= 1 / 3 * sum(x[i, :]) # sum = 3
            y["R", -3, i] >= 1 - 1 / 3 * sum(x[i, :]) # sum = 2
            y["R", -4, i] >= 1 - 1 / 2 * sum(x[i, :]) # sum = 1
            y["R", -5, i] >= 1 - 1 / 1 * sum(x[i, :]) # sum = 0
            # Columns
            y["C", 5, i] <= 1 / 5 * sum(x[:, i]) # sum = 5
            y["C", 4, i] <= 1 / 4 * sum(x[:, i]) # sum = 4
            y["C", 3, i] <= 1 / 3 * sum(x[:, i]) # sum = 3
            y["C", -3, i] >= 1 - 1 / 3 * sum(x[:, i]) # sum = 2
            y["C", -4, i] >= 1 - 1 / 2 * sum(x[:, i]) # sum = 1
            y["C", -5, i] >= 1 - 1 / 1 * sum(x[:, i]) # sum = 0
        end)
    end
    # Solve it
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
    	println("Urban_plan ",solver_name(model)," ", solve_time(model))
    else
    	println("Urban_plan ",solver_name(model)," failed")
    end
    return
end

function example_transp(factory)
    ORIG = ["GARY", "CLEV", "PITT"]
    DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
    supply = [1_400, 2_600, 2_900]
    demand = [900, 1_200, 600, 400, 1_700, 1_100, 1_000]
    cost = [
        39 14 11 14 16 82 8
        27 9 12 9 26 95 17
        24 14 17 13 28 99 20
    ]
    model = Model(factory)
    set_silent(model)
    @variable(model, trans[1:length(ORIG), 1:length(DEST)] >= 0)
    @objective(
        model,
        Min,
        sum(
            cost[i, j] * trans[i, j] for i in 1:length(ORIG),
            j in 1:length(DEST)
        )
    )
    @constraints(
        model,
        begin
            [i in 1:length(ORIG)], sum(trans[i, :]) == supply[i]
            [j in 1:length(DEST)], sum(trans[:, j]) == demand[j]
        end
    )
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
    	println("transp ",solver_name(model)," ", solve_time(model))
    else
    	println("transp ",solver_name(model)," failed")
    end
    return
end

function example_steelT3(factory; verbose = false)
    T = 4
    prod = ["bands", "coils"]
    area = Dict(
        "bands" => ("east", "north"),
        "coils" => ("east", "west", "export"),
    )
    avail = [40, 40, 32, 40]
    rate = Dict("bands" => 200, "coils" => 140)
    inv0 = Dict("bands" => 10, "coils" => 0)
    prodcost = Dict("bands" => 10, "coils" => 11)
    invcost = Dict("bands" => 2.5, "coils" => 3)
    revenue = Dict(
        "bands" => Dict(
            "east" => [25.0, 26.0, 27.0, 27.0],
            "north" => [26.5, 27.5, 28.0, 28.5],
        ),
        "coils" => Dict(
            "east" => [30, 35, 37, 39],
            "west" => [29, 32, 33, 35],
            "export" => [25, 25, 25, 28],
        ),
    )
    market = Dict(
        "bands" => Dict(
            "east" => [2000, 2000, 1500, 2000],
            "north" => [4000, 4000, 2500, 4500],
        ),
        "coils" => Dict(
            "east" => [1000, 800, 1000, 1100],
            "west" => [2000, 1200, 2000, 2300],
            "export" => [1000, 500, 500, 800],
        ),
    )
    # Model
    model = Model(factory)
    set_silent(model)
    # Decision Variables
    @variables(
        model,
        begin
            make[p in prod, t in 1:T] >= 0
            inventory[p in prod, t in 0:T] >= 0
            0 <= sell[p in prod, a in area[p], t in 1:T] <= market[p][a][t]
        end
    )
    @constraints(
        model,
        begin
            [p = prod, a = area[p], t = 1:T], sell[p, a, t] <= market[p][a][t]
            # Total of hours used by all products may not exceed hours available,
            # in each week
            [t in 1:T], sum(1 / rate[p] * make[p, t] for p in prod) <= avail[t]
            # Initial inventory must equal given value
            [p in prod], inventory[p, 0] == inv0[p]
            # Tons produced and taken from inventory must equal tons sold and put
            # into inventory.
            [p in prod, t in 1:T],
            make[p, t] + inventory[p, t-1] ==
            sum(sell[p, a, t] for a in area[p]) + inventory[p, t]
        end
    )
    # Maximize total profit: total revenue less costs for all products in all
    # weeks.
    @objective(
        model,
        Max,
        sum(
            revenue[p][a][t] * sell[p, a, t] - prodcost[p] * make[p, t] -
            invcost[p] * inventory[p, t] for p in prod, a in area[p], t in 1:T
        )
    )
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
    	println("SteelT3 ",solver_name(model)," ", solve_time(model))
    else
    	println("SteelT3 ",solver_name(model)," failed")
    end

    return
end

function example_prod(factory; verbose = false)
    # PRODUCTION SETS AND PARAMETERS
    prd = ["18REG" "24REG" "24PRO"]
    # Members of the product group
    numprd = length(prd)
    pt = [1.194, 1.509, 1.509]
    # Crew-hours to produce 1000 units
    pc = [2304, 2920, 2910]
    # Nominal production cost per 1000, used
    # to compute inventory and shortage costs
    #
    # TIME PERIOD SETS AND PARAMETERS
    firstperiod = 1
    # Index of first production period to be modeled
    lastperiod = 13
    # Index of last production period to be modeled
    numperiods = firstperiod:lastperiod
    # 'planning horizon' := first..last;
    # EMPLOYMENT PARAMETERS
    # Workers per crew
    cs = 18
    # Regular-time hours per shift
    sl = 8
    # Wage per hour for regular-time labor
    rtr = 16.00
    # Wage per hour for overtime labor
    otr = 43.85
    # Crews employed at start of first period
    iw = 8
    # Regular working days in a production period
    dpp = [19.5, 19, 20, 19, 19.5, 19, 19, 20, 19, 20, 20, 18, 18]
    # Maximum crew-hours of overtime in a period
    ol = [96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96]
    # Lower limit on average employment in a period
    cmin = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    # Upper limit on average employment in a period
    cmax = [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
    # Penalty cost of hiring a crew
    hc = [
        7500,
        7500,
        7500,
        7500,
        15000,
        15000,
        15000,
        15000,
        15000,
        15000,
        7500,
        7500,
        7500,
    ]
    # Penalty cost of laying off a crew
    lc = [
        7500,
        7500,
        7500,
        7500,
        15000,
        15000,
        15000,
        15000,
        15000,
        15000,
        7500,
        7500,
        7500,
    ]
    # DEMAND PARAMETERS
    d18REG = [
        63.8,
        76,
        88.4,
        913.8,
        115,
        133.8,
        79.6,
        111,
        121.6,
        470,
        78.4,
        99.4,
        140.4,
        63.8,
    ]
    d24REG = [
        1212,
        306.2,
        319,
        208.4,
        298,
        328.2,
        959.6,
        257.6,
        335.6,
        118,
        284.8,
        970,
        343.8,
        1212,
    ]
    d24PRO = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1102, 0, 0, 0, 0]
    # Requirements (in 1000s) to be met from current production and inventory
    dem = Array[d18REG, d24REG, d24PRO]
    # true if product will be the subject of a special promotion in the period
    pro = Array[
        [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
        [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
    ]
    # INVENTORY AND SHORTAGE PARAMETERS
    # Proportion of non-promoted demand that must be in inventory the previous
    # period
    rir = 0.75
    # Proportion of promoted demand that must be in inventory the previous
    # period
    pir = 0.80
    # Upper limit on number of periods that any product may sit in inventory
    life = 2
    # Inventory cost per 1000 units is cri times nominal production cost
    cri = [0.015, 0.015, 0.015]
    # Shortage cost per 1000 units is crs times nominal production cost
    crs = [1.1, 1.1, 1.1]
    # Inventory at start of first period; age unknown
    iinv = [82, 792.2, 0]
    # Initial inventory still available for allocation at end of period t
    iil = [
        [
            max(0, iinv[p] - sum(dem[p][v] for v in firstperiod:t)) for
            t in numperiods
        ] for p in 1:numprd
    ]
    # Lower limit on inventory at end of period t
    function checkpro(
        product,
        timeperiod,
        production,
        promotionalrate,
        regularrate,
    )
        if production[product][timeperiod+1] == 1
            return promotionalrate
        else
            return regularrate
        end
    end
    minv = [
        [dem[p][t+1] * checkpro(p, t, pro, pir, rir) for t in numperiods]
        for p in 1:numprd
    ]
    # DEFINE MODEL
    prod = Model(factory)
    set_silent(prod)
    # VARIABLES
    # Average number of crews employed in each period
    @variable(prod, Crews[0:lastperiod] >= 0)
    # Crews hired from previous to current period
    @variable(prod, Hire[numperiods] >= 0)
    # Crews laid off from previous to current period
    @variable(prod, Layoff[numperiods] >= 0)
    # Production using regular-time labor, in 1000s
    @variable(prod, Rprd[1:numprd, numperiods] >= 0)
    # Production using overtime labor, in 1000s
    @variable(prod, Oprd[1:numprd, numperiods] >= 0)
    # a numperiods old -- produced in period (t+1)-a --
    # and still in storage at the end of period t
    @variable(prod, Inv[1:numprd, numperiods, 1:life] >= 0)
    # Accumulated unsatisfied demand at the end of period t
    @variable(prod, Short[1:numprd, numperiods] >= 0)
    # CONSTRAINTS
    # Hours needed to accomplish all regular-time production in a period must
    # not exceed hours available on all shifts
    @constraint(
        prod,
        [t = numperiods],
        sum(pt[p] * Rprd[p, t] for p in 1:numprd) <= sl * dpp[t] * Crews[t]
    )
    # Hours needed to accomplish all overtime production in a period must not
    # exceed the specified overtime limit
    @constraint(
        prod,
        [t = numperiods],
        sum(pt[p] * Oprd[p, t] for p in 1:numprd) <= ol[t]
    )
    # Use given initial workforce
    @constraint(prod, Crews[firstperiod-1] == iw)
    # Workforce changes by hiring or layoffs
    @constraint(
        prod,
        [t in numperiods],
        Crews[t] == Crews[t-1] + Hire[t] - Layoff[t]
    )
    # Workforce must remain within specified bounds
    @constraint(prod, [t in numperiods], cmin[t] <= Crews[t])
    @constraint(prod, [t in numperiods], Crews[t] <= cmax[t])
    # 'first demand requirement
    @constraint(
        prod,
        [p in 1:numprd],
        Rprd[p, firstperiod] + Oprd[p, firstperiod] + Short[p, firstperiod] -
        Inv[p, firstperiod, 1] == max(0, dem[p][firstperiod] - iinv[p])
    )
    # Production plus increase in shortage plus decrease in inventory must
    # equal demand
    for t in (firstperiod+1):lastperiod
        @constraint(
            prod,
            [p in 1:numprd],
            Rprd[p, t] + Oprd[p, t] + Short[p, t] - Short[p, t-1] +
            sum(Inv[p, t-1, a] - Inv[p, t, a] for a in 1:life) ==
            max(0, dem[p][t] - iil[p][t-1])
        )
    end
    # Inventory in storage at end of period t must meet specified minimum
    @constraint(
        prod,
        [p in 1:numprd, t in numperiods],
        sum(Inv[p, t, a] + iil[p][t] for a in 1:life) >= minv[p][t]
    )
    # In the vth period (starting from first) no inventory may be more than v
    # numperiods old (initial inventories are handled separately)
    @constraint(
        prod,
        [p in 1:numprd, v in 1:(life-1), a in (v+1):life],
        Inv[p, firstperiod+v-1, a] == 0
    )
    # New inventory cannot exceed production in the most recent period
    @constraint(
        prod,
        [p in 1:numprd, t in numperiods],
        Inv[p, t, 1] <= Rprd[p, t] + Oprd[p, t]
    )
    # Inventory left from period (t+1)-p can only decrease as time goes on
    secondperiod = firstperiod + 1
    @constraint(
        prod,
        [p in 1:numprd, t in 2:lastperiod, a in 2:life],
        Inv[p, t, a] <= Inv[p, t-1, a-1]
    )
    # OBJECTIVE
    # Full regular wages for all crews employed, plus penalties for hiring and
    # layoffs, plus wages for any overtime worked, plus inventory and shortage
    # costs. (All other production costs are assumed to depend on initial
    # inventory and on demands, and so are not included explicitly.)
    @objective(
        prod,
        Min,
        sum(
            rtr * sl * dpp[t] * cs * Crews[t] +
            hc[t] * Hire[t] +
            lc[t] * Layoff[t] +
            sum(
                otr * cs * pt[p] * Oprd[p, t] +
                sum(cri[p] * pc[p] * Inv[p, t, a] for a in 1:life) +
                crs[p] * pc[p] * Short[p, t] for p in 1:numprd
            ) for t in numperiods
        )
    )
    # Obtain solution
    optimize!(prod)
    model = prod
    if termination_status(model) == MOI.OPTIMAL
    	println("prod ",solver_name(model)," ", solve_time(model))
    else
    	println("prod ",solver_name(model)," failed")
    end

    return
end

function example_knapsack(factory; verbose = false)
    profit = [5, 3, 2, 7, 4]
    weight = [2, 8, 4, 2, 5]
    capacity = 10
    model = Model(factory)
    set_silent(model)
    @variable(model, x[1:5], Bin)
    # Objective: maximize profit
    @objective(model, Max, profit' * x)
    # Constraint: can carry all
    @constraint(model, weight' * x <= capacity)
    # Solve problem using MIP solver
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
    	println("knapsack ",solver_name(model)," ", solve_time(model))
    else
    	println("knapsack ",solver_name(model)," failed")
    end
    return
end

function example_factory_schedule(factory)
    # Sets in the problem:
    months, factories = 1:12, [:A, :B]
    # This function takes a matrix and converts it to a JuMP container so we can
    # refer to elements such as `d_max_cap[1, :A]`.
    containerize(A::Matrix) = Containers.DenseAxisArray(A, months, factories)
    # Maximum production capacity in (month, factory) [units/month]:
    d_max_cap = containerize(
        [
            100000 50000
            110000 55000
            120000 60000
            145000 100000
            160000 0
            140000 70000
            155000 60000
            200000 100000
            210000 100000
            197000 100000
            80000 120000
            150000 150000
        ],
    )
    # Minimum production capacity in (month, factory) [units/month]:
    d_min_cap = containerize(
        [
            20000 20000
            20000 20000
            20000 20000
            20000 20000
            20000 0
            20000 20000
            20000 20000
            20000 20000
            20000 20000
            20000 20000
            20000 20000
            20000 20000
        ],
    )
    # Variable cost of production in (month, factory) [$/unit]:
    d_var_cost = containerize([
        10 5
        11 4
        12 3
        9 5
        8 0
        8 6
        5 4
        7 6
        9 8
        10 11
        8 10
        8 12
    ])
    # Fixed cost of production in (month, factory) # [$/month]:
    d_fixed_cost = containerize(
        [
            500 600
            500 600
            500 600
            500 600
            500 0
            500 600
            500 600
            500 600
            500 600
            500 600
            500 600
            500 600
        ],
    )
    # Demand in each month [units/month]:
    d_demand = [
        120_000,
        100_000,
        130_000,
        130_000,
        140_000,
        130_000,
        150_000,
        170_000,
        200_000,
        190_000,
        140_000,
        100_000,
    ]
    # The model!
    model = Model(factory)
    set_silent(model)
    # Decision variables
    @variables(model, begin
        status[m in months, f in factories], Bin
        production[m in months, f in factories], Int
    end)
    # The production cannot be less than minimum capacity.
    @constraint(
        model,
        [m in months, f in factories],
        production[m, f] >= d_min_cap[m, f] * status[m, f],
    )
    # The production cannot be more that maximum capacity.
    @constraint(
        model,
        [m in months, f in factories],
        production[m, f] <= d_max_cap[m, f] * status[m, f],
    )
    # The production must equal demand in a given month.
    @constraint(model, [m in months], sum(production[m, :]) == d_demand[m])
    # Factory B is shut down during month 5, so production and status are both
    # zero.
    fix(status[5, :B], 0.0)
    fix(production[5, :B], 0.0)
    # The objective is to minimize the cost of production across all time
    ##periods.
    @objective(
        model,
        Min,
        sum(
            d_fixed_cost[m, f] * status[m, f] +
            d_var_cost[m, f] * production[m, f] for m in months, f in factories
        )
    )
    # Optimize the problem
    optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
    	println("factory_schedule ",solver_name(model)," ", solve_time(model))
    else
    	println("factory_schedule ",solver_name(model)," failed")
    end
    return
end

for p in [example_factory_schedule, example_knapsack, example_prod, example_steelT3, example_transp, example_urban_plan ]
  for s in [GLPK.Optimizer, with_optimizer(Cbc.Optimizer, logLevel = 0)]
    for trial in 1:5
       p(s)
     end
   end
end