using JuMP
using GLPKMathProgInterface

m = Model(solver = GLPKSolverMIP())
n = 12
M = 1000000

FixCost = [3500, 9000, 10000, 4000, 3000, 9000, 9000, 3000, 4000, 10000, 9000, 3500]
Cap = [300, 250, 100, 180, 275, 300, 200, 220, 280, 250, 230, 180]
Demand = [120, 80, 75, 100, 110, 100, 90, 60, 30, 150, 95, 120]

DeliveryCost = [[100 80 50 50 60 100 120 90 60 70 65 110];
                [120 90 60 70 65 110 140 110 80 80 75 130];
                [140 110 80 80 75 130 160 125 100 100 80 150];
                [160 125 100 100 80 150 190 150 130 M M M];
                [190 150 130 M M M 200 180 150 M M M];
                [200 180 150 M M M 100 80 50 50 60 100];
                [100 80 50 50 60 10 120 90 70 70 65 110];
                [120 90 60 70 65 110 140 110 80 80 75 130];
                [140 110 80 80 75 130 160 125 100 100 80 150];
                [160 125 100 100 80 150 190 150 130 M M M];
                [190 150 130 M M M 200 180 150 M M M];
                [200 180 150 M M M 100 80 50 50 60 100]];

@variable(m, x[1:n,1:n], Bin)

@objective(m, Min, sum(sum(x[i,j]*DeliveryCost[i,j] + FixCost[j]  for j = 1:n) for i = 1:n))

@constraint(m, [i=1:n], sum(x[i,j] == 1 for j = 1:n))
@constraint(m, [j=1:n], sum(x[i,j] <= 4 for i = 1:n))
@constraint(m, [j=1:n], sum(x[i,j]*Demand[i] for i = 1:n) <= Cap[j])
@constraint(m, [i=1:n], x[i,5]+x[i,6] == 1)

status = solve(m)

println("Objective = ", getobjectivevalue(m))
println("X = ", getvalue(x))