CMCLass: 派送問題(2)

如果派送問題(1)改成每個軍事基地最多可容納2個舊空軍基地，那麼怎麼算呢?

 1  #!/usr/bin/env python
 2  # -*- coding: utf-8 -*-
 3
 4  def caldest(From, To):
 5      from math import sqrt
 6      dest = []
 7
 8      for start in From:
 9          subdest = []
10          for end in To:
11              subdest.append(round(sqrt((start[0] - end[0])**2 + (start[1] - end[1])**2), 0))
12          dest.append(subdest)
13
14      return dest
15
16  def optima(coef):
17      from pulp import *
18      (Dest, var, obj) = (coef, [], 0)
19      prob = LpProblem("dest", LpMinimize)
20
21      for i in xrange(len(Dest)):
22          var.append([])
23          for j in xrange(len(Dest[0])):
24              (str_i, str_j) = ("%03d" % (i + 1), "%03d" % (j + 1))
25              var[i].append(LpVariable("var"+str_i+'_'+str_j, 0, 1, LpInteger))
26              obj += Dest[i][j] * var[i][j]
27
28      prob += obj, 'OBJ.'
29
30      for i in xrange(len(Dest)):
31          st = 0
32          for j in xrange(len(Dest[0])):
33              st += var[i][j]
34          prob += st == 1
35
36      for j in xrange(len(Dest[0])):
37          st = 0
38          for i in xrange(len(Dest)):
39              st += var[i][j]
40          prob += st <= 2
41
42      #prob.writeLP("dest.lp")
43      prob.solve()
44
45      print "Status:", LpStatus[prob.status]
46
47      for v in prob.variables():
48          if v.varValue >= 1:
49              print v.name[3:6] + ' => ' + v.name[7:10]
50
51      print "objective=", value(prob.objective)
52
53  if __name__ == '__main__':
54      A = [[1130, 863],
55          [1705, 1283],
56          [1326, 1736],
57          [975, 825],
58          [1286, 807],
59          [909, 1143],
60          [1579, 1608],
61          [1162, 1118],
62          [876, 1573],
63          [1198, 1748],
64          [801, 1134],
65          [956, 1226]]
66
67      B = [[1031, 1206],
68          [1046, 1000],
69          [1803, 1809],
70          [1588, 1821],
71          [1300, 1482],
72          [1939, 1228],
73          [1147, 1703],
74          [1319, 1254],
75          [1300, 1182],
76          [1817, 1269],
77          [1945, 1605],
78          [1349, 1488]]
79
80      dest = caldest(A, B)
81      optima(dest)

在第40行的地方改成 <= 2 即可