In [14]:
#容量制約付き施設配置問題
from gurobipy import *
def make_data():
    I,d = multidict({1:80, 2:270, 3:250, 4:160, 5:180})            # demand
    J,M,f = multidict({1:[500,1000], 2:[500,1000], 3:[500,1000]}) # capacity, fixed costs
    c = {(1,1):4,  (1,2):6,  (1,3):9,    # transportation costs
         (2,1):5,  (2,2):4,  (2,3):7,
         (3,1):6,  (3,2):3,  (3,3):4,
         (4,1):8,  (4,2):5,  (4,3):3,
         (5,1):10, (5,2):8,  (5,3):4,
         }
    return I,J,d,M,f,c

In [15]:
#multidictを三つ一気にやると、keyの先のリストから一つずつ辞書の先にくっつく
make_data()


Out[15]:
([1, 2, 3, 4, 5],
 [1, 2, 3],
 {1: 80, 2: 270, 3: 250, 4: 160, 5: 180},
 {1: 500, 2: 500, 3: 500},
 {1: 1000, 2: 1000, 3: 1000},
 {(1, 1): 4,
  (1, 2): 6,
  (1, 3): 9,
  (2, 1): 5,
  (2, 2): 4,
  (2, 3): 7,
  (3, 1): 6,
  (3, 2): 3,
  (3, 3): 4,
  (4, 1): 8,
  (4, 2): 5,
  (4, 3): 3,
  (5, 1): 10,
  (5, 2): 8,
  (5, 3): 4})

In [39]:
#変数の名前はつけないと結果の解釈ができない
def flp(I,J,d,M,f,c):
    model = Model("flp")
    x, y = {} ,{}
    for j in J:
        y[j] = model.addVar(vtype = "B", name = "facility(%s)" %j)
        for i in I:
            x[i,j] = model.addVar(vtype = "C", name = "transport(%s, %s)" %(i,j))
    model.update()
    
    for i in I:
        model.addConstr(quicksum(x[i,j] for j in J) == d[i])
    for j in J:
        model.addConstr(quicksum(x[i,j] for i in I) <= M[j]* y[j])
    for (i,j) in x:
        model.addConstr(x[i,j] <= d[i] * y[j])
        
    model.setObjective(quicksum(f[j]*y[j] for j in J) + quicksum(c[i,j]* x[i,j] for i in I for j in J))
    model.__data = x,y
    
    return model

In [40]:
if __name__ == "__main__":
    I,J,d,c,f,M = make_data()
    model = flp(I,J,d,c,f,M)
    model.optimize()


Optimize a model with 23 rows, 18 columns and 63 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 5e+02]
  Objective range [3e+00, 1e+03]
  Bounds range    [1e+00, 1e+00]
  RHS range       [8e+01, 3e+02]
Found heuristic solution: objective 7990
Presolve time: 0.00s
Presolved: 23 rows, 18 columns, 63 nonzeros
Variable types: 15 continuous, 3 integer (3 binary)

Root relaxation: objective 5.610000e+03, 7 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0    5610.0000000 5610.00000  0.00%     -    0s

Explored 0 nodes (7 simplex iterations) in 0.01 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 5.610000000000e+03, best bound 5.610000000000e+03, gap 0.0%

In [41]:
for v in model.getVars():
    print v.VarName ,v.X


facility(1) -0.0
transport(1, 1) 0.0
transport(2, 1) 0.0
transport(3, 1) 0.0
transport(4, 1) 0.0
transport(5, 1) 0.0
facility(2) 1.0
transport(1, 2) 80.0
transport(2, 2) 270.0
transport(3, 2) 150.0
transport(4, 2) 0.0
transport(5, 2) 0.0
facility(3) 1.0
transport(1, 3) 0.0
transport(2, 3) 0.0
transport(3, 3) 100.0
transport(4, 3) 160.0
transport(5, 3) 180.0

In [53]:
#k-median problem
def kmedian(I, J, c, k):
    model = Model("k-median")
    x ,y = {} ,{}
    for j in J:
        y[j] = model.addVar(vtype= "B", name = "facility(%s)" %j)
        
        for i in I:
            x[i,j] = model.addVar(vtype="B" ,name = "satisfaction(%s, %s)" %(i,j))
    model.update()
    
    for i in I:
        model.addConstr(quicksum(x[i,j] for j in J) == 1)
        for j in J:
            model.addConstr(x[i, j] <= y[j])
    
    model.addConstr(quicksum(y[j] for j in J) == k)
    model.setObjective(quicksum(c[i,j] *x[i,j] for i in I for j in J))
    
    model.__data = x, y
    return model

In [54]:
def make_data2():
    I = [1,2,3,4,5]
    J = [1,2,3]
    c = c = {(1,1):4,  (1,2):6,  (1,3):9, 
         (2,1):5,  (2,2):4,  (2,3):7,
         (3,1):6,  (3,2):3,  (3,3):4,
         (4,1):8,  (4,2):5,  (4,3):3,
         (5,1):10, (5,2):8,  (5,3):4,}
    k = 2
    return I,J,c,k

In [55]:
if __name__ == "__main__":
    I,J,c,k = make_data2()
    model = kmedian( I,J,c,k)
    model.optimize()


Optimize a model with 21 rows, 18 columns and 48 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [3e+00, 1e+01]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 2e+00]
Found heuristic solution: objective 27
Presolve time: 0.00s
Presolved: 21 rows, 18 columns, 48 nonzeros
Variable types: 0 continuous, 18 integer (18 binary)

Root relaxation: objective 2.000000e+01, 11 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0      20.0000000   20.00000  0.00%     -    0s

Explored 0 nodes (11 simplex iterations) in 0.01 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

In [56]:
for v in model.getVars():
    print v.VarName, v.X
    
#1、2さんは施設1に
#3、4、5さんは施設3に
#施設2は開設されなかった


facility(1) 1.0
satisfaction(1, 1) 1.0
satisfaction(2, 1) 1.0
satisfaction(3, 1) -0.0
satisfaction(4, 1) -0.0
satisfaction(5, 1) -0.0
facility(2) 0.0
satisfaction(1, 2) -0.0
satisfaction(2, 2) 0.0
satisfaction(3, 2) 0.0
satisfaction(4, 2) -0.0
satisfaction(5, 2) -0.0
facility(3) 1.0
satisfaction(1, 3) -0.0
satisfaction(2, 3) -0.0
satisfaction(3, 3) 1.0
satisfaction(4, 3) 1.0
satisfaction(5, 3) 1.0

In [57]:
#k-center problem
def kcenter(I, J, c, k):
    model = Model("k-center")
    z = model.addVar(vtype = "C", name = "max_distance" )
    x, y = {}, {}
    
    for j in J:
        y[j] = model.addVar(vtype= "B", name = "facility(%s)" %j)
        
        for i in I:
            x[i,j] = model.addVar(vtype="B" ,name = "satisfaction(%s, %s)" %(i,j))
    model.update()
    
    for i in I:
        model.addConstr(quicksum(x[i,j] for j in J) == 1)
        for j in J:
            model.addConstr(x[i, j] <= y[j])
            model.addConstr(c[i,j] * x[i,j] <= z)
    
    model.addConstr(quicksum(y[j] for j in J) == k)
    model.setObjective(z)
    model.__data = x, y
    return model

In [58]:
if __name__ == "__main__":
    I,J,c,k = make_data2()
    model = kcenter( I,J,c,k)
    model.optimize()


Optimize a model with 36 rows, 19 columns and 78 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+01]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 2e+00]
Found heuristic solution: objective 8
Presolve time: 0.00s
Presolved: 36 rows, 19 columns, 78 nonzeros
Variable types: 1 continuous, 18 integer (18 binary)

Root relaxation: objective 4.000000e+00, 9 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    4.00000    0    8    8.00000    4.00000  50.0%     -    0s
H    0     0                       6.0000000    4.00000  33.3%     -    0s
H    0     0                       5.0000000    4.00000  20.0%     -    0s

Explored 0 nodes (10 simplex iterations) in 0.01 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 5.000000000000e+00, best bound 5.000000000000e+00, gap 0.0%

In [59]:
for v in model.getVars():
    print v.VarName, v.X
    
#この場合はk-median problemと解は同じ


max_distance 5.0
facility(1) 1.0
satisfaction(1, 1) 1.0
satisfaction(2, 1) 1.0
satisfaction(3, 1) -0.0
satisfaction(4, 1) -0.0
satisfaction(5, 1) -0.0
facility(2) -0.0
satisfaction(1, 2) -0.0
satisfaction(2, 2) -0.0
satisfaction(3, 2) -0.0
satisfaction(4, 2) -0.0
satisfaction(5, 2) -0.0
facility(3) 1.0
satisfaction(1, 3) -0.0
satisfaction(2, 3) -0.0
satisfaction(3, 3) 1.0
satisfaction(4, 3) 1.0
satisfaction(5, 3) 1.0

In [ ]: