数学规划-指派问题_一次指派问题与二次指派问题

问题:

某公司指派n个员工到n个城市工作(每个城市单独一人),希望使所花费的总电话费用
尽可能少.n个员工两两之间每个月通话的时间表示在下面的矩阵的上三角形部分(假
设通话的时间矩阵是对称的,没有必要写出下三角形部分),n个城市两两之间通话费率
表示在下面的矩阵的下三角形部分(同样道理,假设通话的费率矩阵是对称的,没有必
要写出上三角形部分)．试求解该二次指派问题.(如果你的软件解不了这么大规模的
问题.那就只考虑最前面的若干员工和城市.)

[0 , 5 , 3 , 7 ,  9 , 3 , 9 , 2 , 9 , 0 ]
[7 , 0 , 7 , 8 ,  3 , 2 , 3 , 3 , 5 , 7 ]
[4 , 8 , 0 , 9 ,  3 , 5 , 3 , 3 , 9 , 3 ]
[6 , 2 , 10, 0 ,  8 , 4 , 1 , 8 , 0 , 4 ]
[8 , 6 , 4 , 6 ,  0 , 8 , 8 , 7 , 5 , 9 ]
[8 , 5 , 4 , 6 ,  6 , 0 , 4 , 8 , 0 , 3 ]
[8 , 6 , 7 , 9 ,  4 , 3 , 0 , 7 , 9 , 5 ]
[6 , 8 , 2 , 3 ,  8 , 8 , 6 , 0 , 5 , 5 ]
[6 , 3 , 6 , 2 ,  8 , 3 , 7 , 8 , 0 , 5 ]
[5 , 6 , 7 , 6 ,  6 , 2 , 8 , 8 , 9 , 0 ]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

设0-1变量$x_{ij}$若为1则代表将第i个员工派向第j个城市，若为0则代表未指派$tim{e}_{im}$表示第i个和第m个员工之间的通话和时间，$pric{e}_{jn}$表示j和n城之间的通话费率，那么即可建立模型如下
M i n q q u a d y = ∑ i , j , m , n = 1 i , j , m , n = 10 x i j ∗ x m n ∗ t i m e i m ∗ p r i c e j n 2 S . T . { ∑ i = 1 i = 10 x i j = = 1 ∑ j = 1 j = 10 x i j = = 1 x i j = { 1 0 Min \\qquad y =\frac{\sum_{i,j,m,n=1}^{i,j,m,n=10} x_{ij}*x_{mn}*time_{im}*price_{jn}}{2}\\ S.T.\left\{

$$\begin{aligned} \sum_{i=1}^{i=10}x_{ij}==1\\ \sum_{j=1}^{j=10}x_{ij}==1\\ x_{ij}=\left\{ \begin{aligned} 1 \\ 0 \\ \end{aligned}$$ \right. \end{aligned} \right. Minqquady=2∑i,j,m,n=1i,j,m,n=10​xij​∗xmn​∗timeim​∗pricejn​​S.T.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​i=1∑i=10​xij​==1j=1∑j=10​xij​==1xij​={10​​

import numpy as np
import gurobipy as grb
from gurobipy import *

rec = np.array([[0 , 5 , 3 , 7 ,  9 , 3 , 9 , 2 , 9 , 0 ]
               ,[7 , 0 , 7 , 8 ,  3 , 2 , 3 , 3 , 5 , 7 ]
               ,[4 , 8 , 0 , 9 ,  3 , 5 , 3 , 3 , 9 , 3 ]
               ,[6 , 2 , 10, 0 ,  8 , 4 , 1 , 8 , 0 , 4 ]
               ,[8 , 6 , 4 , 6 ,  0 , 8 , 8 , 7 , 5 , 9 ]
               ,[8 , 5 , 4 , 6 ,  6 , 0 , 4 , 8 , 0 , 3 ]
               ,[8 , 6 , 7 , 9 ,  4 , 3 , 0 , 7 , 9 , 5 ]
               ,[6 , 8 , 2 , 3 ,  8 , 8 , 6 , 0 , 5 , 5 ]
               ,[6 , 3 , 6 , 2 ,  8 , 3 , 7 , 8 , 0 , 5 ]
               ,[5 , 6 , 7 , 6 ,  6 , 2 , 8 , 8 , 9 , 0 ]])
triu = np.triu(rec)#取出矩阵的上三角部分
tril = np.tril(rec)#取出矩阵的下三角部分
time = triu+triu.T
price = tril+tril.T
print("各员工通话时间为:\n",time,"\n","各城市通话单价为:\n",price)

model = grb.Model("Company")
x = model.addVars(10,10,vtype=grb.GRB.BINARY,name = "是否指派")
objective = grb.quicksum(x[i,j] * x[m,n] * time[i,m] * price[j,n] for i in range(10) for j in range(10) for m in range(10) for n in range(10))
#设置目标函数
model.setObjective(objective/2,grb.GRB.MINIMIZE)
#添加约束
model.addConstrs((x.sum(i,'*')==1 for i in range(10)),"行")
model.addConstrs((x.sum('*',i)==1 for i in range(10)),"列")
# 求解
model.optimize()
print('目标函数值是：', model.objVal)
if model.status == GRB.OPTIMAL:
    model.printAttr('X')
    model.printAttr('Slack')
    print("指派矩阵为:\n",np.array(model.X).reshape(10,10))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

各员工通话时间为:
 [[0 5 3 7 9 3 9 2 9 0]
 [5 0 7 8 3 2 3 3 5 7]
 [3 7 0 9 3 5 3 3 9 3]
 [7 8 9 0 8 4 1 8 0 4]
 [9 3 3 8 0 8 8 7 5 9]
 [3 2 5 4 8 0 4 8 0 3]
 [9 3 3 1 8 4 0 7 9 5]
 [2 3 3 8 7 8 7 0 5 5]
 [9 5 9 0 5 0 9 5 0 5]
 [0 7 3 4 9 3 5 5 5 0]] 
 各城市通话单价为:
 [[ 0  7  4  6  8  8  8  6  6  5]
 [ 7  0  8  2  6  5  6  8  3  6]
 [ 4  8  0 10  4  4  7  2  6  7]
 [ 6  2 10  0  6  6  9  3  2  6]
 [ 8  6  4  6  0  6  4  8  8  6]
 [ 8  5  4  6  6  0  3  8  3  2]
 [ 8  6  7  9  4  3  0  6  7  8]
 [ 6  8  2  3  8  8  6  0  8  8]
 [ 6  3  6  2  8  3  7  8  0  9]
 [ 5  6  7  6  6  2  8  8  9  0]]
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 20 rows, 100 columns and 200 nonzeros
Model fingerprint: 0xc0a57822
Model has 3780 quadratic objective terms
Variable types: 0 continuous, 100 integer (100 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [0e+00, 0e+00]
  QObjective range [4e+00, 2e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 1397.0000000
Presolve time: 0.01s
Presolved: 20 rows, 100 columns, 200 nonzeros
Presolved model has 3880 quadratic objective terms
Variable types: 0 continuous, 100 integer (100 binary)

Root relaxation: objective -3.181574e+03, 120 iterations, 0.02 seconds

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

H    0     0                    1329.0000000    0.00000   100%     -    0s
     0     0   -0.00000    0  100 1329.00000    0.00000   100%     -    0s
H    0     0                    1306.0000000    0.00000   100%     -    0s
H    0     0                    1267.0000000    0.00000   100%     -    0s
     0     0   -0.00000    0  100 1267.00000    0.00000   100%     -    0s
H    0     0                    1264.0000000    0.00000   100%     -    0s
     0     2   -0.00000    0  100 1264.00000    0.00000   100%     -    0s
H   37    40                    1261.0000000    0.00000   100%   3.3    0s
H   71    88                    1252.0000000    0.00000   100%   3.1    0s
H 1173   816                    1249.0000000    2.00000   100%   7.3    1s
*15853  7601              56    1239.0000000  230.86615  81.4%  10.7    4s
 18849  8805  536.85272   46   27 1239.00000  248.95813  79.9%  10.6    5s
 30142 13049  321.91415   39  150 1239.00000  321.91415  74.0%  10.2   10s
H30176 12420                    1236.0000000  321.91415  74.0%  10.2   11s
*31767 12255              71    1235.0000000  425.65511  65.5%  11.1   12s
H32185 11796                    1231.0000000  434.92542  64.7%  11.3   12s
H35370 11959                    1225.0000000  626.06833  48.9%  12.3   13s
 37995 12792 1082.77385   59   43 1225.00000  646.16346  47.3%  12.9   15s
*49553 13460              61    1210.0000000  804.21558  33.5%  14.3   18s
*51315 12861              57    1194.0000000  814.04601  31.8%  14.3   18s
*54842 11950              58    1181.0000000  825.96584  30.1%  14.5   19s
 57575 11873 1068.54293   56   46 1181.00000  835.93457  29.2%  14.6   20s
 81670 17253 1178.88719   59   31 1181.00000  888.77588  24.7%  15.2   25s
 107671 20010 1113.69678   61   38 1181.00000  991.67959  16.0%  15.2   30s
H109601 16775                    1142.0000000  994.64106  12.9%  15.2   30s
 132285  8001     cutoff   56      1142.00000 1054.99461  7.62%  15.1   35s

Cutting planes:
  Gomory: 4
  Implied bound: 15
  MIR: 15
  StrongCG: 1
  Flow cover: 103
  Zero half: 1

Explored 143305 nodes (2134002 simplex iterations) in 36.46 seconds
Thread count was 8 (of 8 available processors)

Solution count 10: 1142 1181 1194 ... 1249

Optimal solution found (tolerance 1.00e-04)
Best objective 1.142000000000e+03, best bound 1.142000000000e+03, gap 0.0000%
目标函数值是： 1142.0

    Variable            X 
-------------------------
   /&>[0,8]            1 
   /&>[1,0]            1 
   /&>[2,7]            1 
   /&>[3,2]            1 
   /&>[4,5]            1 
   /&>[5,6]            1 
   /&>[6,1]            1 
   /&>[7,4]            1 
   /&>[8,3]            1 
   /&>[9,9]            1 

  Constraint        Slack 
-------------------------
指派矩阵为:
 [[0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
 [1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115