TSP问题的优化求解

TSP旅行商问题

1.原题

现有一名外国游客从广州入境，他想在144小时以内游玩尽可能多的城市，同时要求综合游玩体验最好，请你规划他的游玩路线。需要结合游客的要求给出具体的游玩路线，包括总花费时间，门票和交通的总费用以及可以游玩的景点数量。他的要求有： （不考虑money）。

为便于测试，这里简化题目：现在有一名游客想在144小时内尽可能多游玩一些城市，给出去往每个城市的时间花费矩阵，进行优化求解得到路线图。

2.模型

优化的目标在除了起点之外的49个城市中找到

Max∑i=0n−1∑j=0n−1xij Max\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} x_{ij}

其中xijx_{ij}为决策变量，表示是否决定从城市ii到城市jj，如果选择去则xij=1x_{ij}=1，否则xij=0x_{ij}=0。

约束条件∑i=0n−1∑j=0n−1Gijxij≤max_time(总时间花费)\text{约束条件} \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} G_{ij} x_{ij} \leq \text{max\_time} \quad \text{(总时间花费)}

GG是输入的时间矩阵，GijG_{ij}表示从城市ii到城市jj的时间花费，max_timemax\_time表示最大时间限制。这里的时间是144−在广州的2.5小时144 -在广州的2.5小时

∑j=0n−1,j≠ixij=∑j=0n−1,j≠ixji∀i∈{0,…,n−1}(流量平衡，进出相等) \sum_{j=0}^{n-1, j \neq i} x_{ij} = \sum_{j=0}^{n-1, j \neq i} x_{ji} \quad \forall i \in \{0, \ldots, n-1\} \quad \text{(流量平衡，进出相等)}

∑j=0n−1,j≠ixij≤1∀i∈{0,…,n−1}(每个城市最多有一个出边) \sum_{j=0}^{n-1, j \neq i} x_{ij} \leq 1 \quad \forall i \in \{0, \ldots, n-1\} \quad \text{(每个城市最多有一个出边)}

∑j=0n−1,j≠ixji≤1∀i∈{0,…,n−1}(每个城市最多有一个入边) \sum_{j=0}^{n-1, j \neq i} x_{ji} \leq 1 \quad \forall i \in \{0, \ldots, n-1\} \quad \text{(每个城市最多有一个入边)}

∑j=0n−1,j≠start_cityxstart_city,j=1(将起点城市出边设置为1表示从这里出发) \sum_{j=0}^{n-1, j \neq \text{start\_city}} x_{\text{start\_city}, j} = 1 \quad \text{(将起点城市出边设置为1表示从这里出发)}

ui−uj+nxij≤n−1∀i,j∈{1,…,n−1},i≠j(MTZ约束，确保不被重复访问) u_i - u_j + n x_{ij} \leq n - 1 \quad \forall i, j \in \{1, \ldots, n-1\}, i \neq j \quad \text{(MTZ约束，确保不被重复访问)}

这里的n是城市的数量，n=50n=50。

对于MTZ约束，这里的uiu_i是一个辅助变量，表示城市ii的访问顺序。如果xij=1x_{ij}=1，那么ui−uj≤−1u_i - u_j \leq -1，即ui≤uj−1<uju_i \leq u_j-1 < u_j.这样就可以保证经过的u是增长的，即保证了城市的访问顺序。

假设存在被重复访问的节点v1v2v3v4,其对应的u为u1u2u3u1,有u1<u2<u3<u1v_1 v_2 v_3 v_4,其对应的u为u_1 u_2 u_3 u_1 , 有u_1<u_2<u_3<u_1产生矛盾，所以会避免这种情况发生。

xii=0∀i∈{0,…,n−1}(防止自环，即防止原地绕圈) x_{ii} = 0 \quad \forall i \in \{0, \ldots, n-1\} \quad \text{(防止自环，即防止原地绕圈)}

变量xij∈{0,1}∀i,j∈{0,…,n−1}\text{变量} \quad x_{ij} \in \{0, 1\} \quad \forall i, j \in \{0, \ldots, n-1\}

ui∈R∀i∈{0,…,n−1} u_i \in \mathbb{R} \quad \forall i \in \{0, \ldots, n-1\}

3.1 Python Gurobipy求解

这里生成300阶的随机数矩阵作为输入。

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
import gurobipy as gp
from gurobipy import GRB


def solve_max_cities_with_time_constraint(Gmatrix, max_time, start_city):
    n = len(Gmatrix)
    '''
    输入参数：时间矩阵 Gmatrix，最大时间 max_time，起点城市 start_city
    输出参数：最佳路径 path，花费时间 total_time，访问的城市数量 cities_visited
    '''
    model = gp.Model("Max_Cities_with_Time_Constraint")

    # 决策变量：x[i, j] ,是否选择从城市 i 到城市 j 的路径
    x = model.addVars(n, n, vtype=GRB.BINARY, name="x")

    # 辅助变量：u[i] 用于避免子环
    u = model.addVars(n, vtype=GRB.CONTINUOUS, name="u")

    # 最大化访问的城市数量
    model.setObjective(gp.quicksum(x[i, j] for i in range(n) for j in range(n)), GRB.MAXIMIZE)

    # 总时间花费 <= max_time
    model.addConstr(gp.quicksum(Gmatrix[i][j] * x[i, j] for i in range(n) for j in range(n)) <= max_time, "time_limit")

    # 添加流量约束
    for i in range(0,n-1):
        model.addConstr(gp.quicksum(x[i, j] for j in range(n) if j != i) == gp.quicksum(x[j, i] for j in range(n) if j != i), f"flow_balance_{i}")
        model.addConstr(gp.quicksum(x[i, j] for j in range(n) if j != i) <= 1, f"out_{i}")
        model.addConstr(gp.quicksum(x[j, i] for j in range(n) if j != i) <= 1, f"in_{i}")

    # 添加起点约束：确保从指定的起点城市出发
    model.addConstr(gp.quicksum(x[start_city, j] for j in range(n) if j != start_city) == 1, "start_city_out")
    #model.addConstr(gp.quicksum(x[i, start_city] for i in range(n) if i != start_city) == 0, "start_city_in")

    # 添加子环消除约束（MTZ约束）
    for i in range(1, n):
        for j in range(1, n):
            if i != j:
                model.addConstr(u[i] - u[j] + n * x[i, j] <= n - 1, f"mtz_{i}_{j}")

    # 防止自环
    for i in range(n):
        model.addConstr(x[i, i] == 0, f"no_self_loop_{i}")

    # 优化模型
    model.optimize()

    # 检查是否找到最优解
    if model.status == GRB.OPTIMAL:
        path = []
        for i in range(n):
            for j in range(n):
                if x[i, j].X > 0.5:  # 判断变量值是否为1
                    path.append((i, j))
        total_time = gp.quicksum(Gmatrix[i][j] * x[i, j].X for i in range(n) for j in range(n)).getValue()
        cities_visited = len(set([i for i, j in path] + [j for i, j in path]))
        return path, total_time, cities_visited
    else:
        print("没有找到最优解")
        return None, None, 0

# 用时间矩阵求解最大城市数目
import numpy as np
Gmatrix = np.random.randint(1, 1000, (300, 300))
max_time = 144 - 2.5 #144减去在广州的2.5小时
path, total_time, cities_visited = solve_max_cities_with_time_constraint(Gmatrix, max_time,4)
print("最佳路径:", path)
print("花费时间:", total_time+2.5)
print("访问的城市数量:", cities_visited)

300阶的求解结果，用时645s：

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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26120.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 90301 rows, 90300 columns and 715509 nonzeros
Model fingerprint: 0x96df0c90
Variable types: 300 continuous, 90000 integer (90000 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+03]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 3e+02]
Presolve removed 318 rows and 77810 columns
Presolve time: 3.95s
Presolved: 89983 rows, 12490 columns, 249326 nonzeros
Variable types: 299 continuous, 12191 integer (12191 binary)

Deterministic concurrent LP optimizer: primal and dual simplex (primal and dual model)
Showing primal log only...

Root relaxation presolved: 12490 rows, 102174 columns, 261517 nonzeros


Root simplex log...

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    4.3750000e+02   0.000000e+00   2.805988e+04      5s
Concurrent spin time: 0.09s

Solved with dual simplex

Root relaxation: objective 6.434205e+01, 1347 iterations, 0.86 seconds (0.28 work units)

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

     0     0   64.34205    0   34          -   64.34205      -     -    5s
     0     0   64.03846    0   31          -   64.03846      -     -   11s
     0     0   63.92804    0   31          -   63.92804      -     -   16s
     0     0   63.92804    0   41          -   63.92804      -     -   17s
     0     0   63.92804    0   41          -   63.92804      -     -   21s
H    0     0                       7.0000000   63.92804   813%     -   21s
     0     0   63.92804    0   33    7.00000   63.92804   813%     -   21s
H    0     0                       9.0000000   63.92804   610%     -   44s
     0     0   63.92804    0   33    9.00000   63.92804   610%     -   44s
     0     0   63.92804    0   39    9.00000   63.92804   610%     -   45s
     0     0   63.92804    0   39    9.00000   63.92804   610%     -   53s
     0     0   63.92804    0   39    9.00000   63.92804   610%     -   53s
     0     0   63.88293    0   49    9.00000   63.88293   610%     -   57s
     0     0   63.88293    0   49    9.00000   63.88293   610%     -   58s
     0     0   63.88293    0   49    9.00000   63.88293   610%     -   62s
     0     0   63.88293    0   49    9.00000   63.88293   610%     -   62s
     0     0   63.88293    0   33    9.00000   63.88293   610%     -   66s
     0     0   63.88293    0   39    9.00000   63.88293   610%     -   67s
H    0     0                      19.0000000   63.88293   236%     -   80s
     0     0   63.88293    0   39   19.00000   63.88293   236%     -   80s
     0     0   63.88293    0   39   19.00000   63.88293   236%     -   81s
     0     0   63.88293    0   39   19.00000   63.88293   236%     -   85s
     0     0   63.88293    0   41   19.00000   63.88293   236%     -   86s
     0     0   63.88293    0   41   19.00000   63.88293   236%     -   89s
     0     0   63.88293    0   41   19.00000   63.88293   236%     -   90s
     0     0   63.84848    0   59   19.00000   63.84848   236%     -   94s
     0     0   63.84848    0   59   19.00000   63.84848   236%     -   95s
     0     0   63.84848    0   59   19.00000   63.84848   236%     -   98s
     0     0   63.84848    0   59   19.00000   63.84848   236%     -   99s
     0     0   63.84848    0   60   19.00000   63.84848   236%     -  102s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  103s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  105s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  106s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  108s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  109s
     0     0   63.84848    0   61   19.00000   63.84848   236%     -  110s
     0     0   63.81818    0   34   19.00000   63.81818   236%     -  118s
     0     0   63.76923    0   41   19.00000   63.76923   236%     -  125s
     0     0   63.66450    0   39   19.00000   63.66450   235%     -  135s
     0     0   63.66450    0   39   19.00000   63.66450   235%     -  135s
     0     0   63.66450    0   39   19.00000   63.66450   235%     -  136s
     0     0   63.66450    0   39   19.00000   63.66450   235%     -  140s
     0     0   63.66450    0   49   19.00000   63.66450   235%     -  140s
     0     0   63.66450    0   49   19.00000   63.66450   235%     -  140s
     0     0   63.66450    0   49   19.00000   63.66450   235%     -  141s
     0     0   63.66450    0   60   19.00000   63.66450   235%     -  144s
     0     0   63.66450    0   61   19.00000   63.66450   235%     -  144s
     0     0   63.66450    0   61   19.00000   63.66450   235%     -  145s
     0     0   63.66450    0   54   19.00000   63.66450   235%     -  145s
H    0     0                      21.0000000   63.66450   203%     -  150s
     0     0   63.66450    0   54   21.00000   63.66450   203%     -  150s
     0     0   63.66450    0   54   21.00000   63.66450   203%     -  151s
     0     2   63.00000    0   54   21.00000   63.00000   200%     -  163s
     3     8   63.00000    2   51   21.00000   63.00000   200%  76.3  167s
    11    16   62.89189    4   44   21.00000   63.00000   200%  60.7  171s
    23    28   62.85714    6   46   21.00000   63.00000   200%  55.7  175s
    35    40   62.78824    8   40   21.00000   63.00000   200%  60.7  180s
    47    52   62.33333   11   25   21.00000   63.00000   200%  54.0  249s
H   49    52                      48.0000000   63.00000  31.3%  52.8  249s
    51    56   62.31606   12   59   48.00000   63.00000  31.3%  52.1  252s
H   52    56                      50.0000000   63.00000  26.0%  51.1  252s
    63    72   62.30566   15   71   50.00000   63.00000  26.0%  46.6  255s
    92    96   62.29142   22   72   50.00000   63.00000  26.0%  38.7  260s
   136   121   62.27479   29   55   50.00000   63.00000  26.0%  35.9  265s
   164   159   62.24837   33   57   50.00000   63.00000  26.0%  52.0  271s
   194   186   62.23794   38   54   50.00000   63.00000  26.0%  50.6  275s
   253   245   62.21036   48   53   50.00000   63.00000  26.0%  43.0  281s
   285   275   62.18348   55   57   50.00000   63.00000  26.0%  40.1  287s
   303   298   62.16999   60   41   50.00000   63.00000  26.0%  39.7  291s
   327   299   61.82886   34   61   50.00000   63.00000  26.0%  39.0  313s
   329   300   61.74333   58   34   50.00000   63.00000  26.0%  38.8  319s
   330   301   62.18166   55   39   50.00000   63.00000  26.0%  38.7  321s
   332   302   60.30097   17   54   50.00000   63.00000  26.0%  38.4  325s
   335   304   63.00000   60   44   50.00000   63.00000  26.0%  38.1  333s
   337   306   61.80000    3   70   50.00000   63.00000  26.0%  37.9  339s
   339   307   61.66667    9   45   50.00000   63.00000  26.0%  37.6  340s
   340   308   59.91273   25   45   50.00000   63.00000  26.0%  37.5  345s
   342   309   63.00000   25   66   50.00000   63.00000  26.0%  37.3  352s
   344   310   62.14956   10   85   50.00000   63.00000  26.0%  37.1  358s
   345   311   56.91556   65   67   50.00000   63.00000  26.0%  37.0  360s
   347   312   59.86378   32   93   50.00000   63.00000  26.0%  36.8  366s
   349   314   60.26698   19   74   50.00000   63.00000  26.0%  36.6  372s
   350   314   61.25759   27   86   50.00000   63.00000  26.0%  36.5  375s
   353   316   61.99417   17   83   50.00000   63.00000  26.0%  36.2  380s
   356   318   61.14295   12   93   50.00000   63.00000  26.0%  35.8  400s
   359   320   61.85048   27  103   50.00000   63.00000  26.0%  35.5  405s
   360   321   61.14295   12  103   50.00000   62.99272  26.0%  35.5  424s
   361   325   62.99272   10   58   50.00000   62.99272  26.0%  45.6  426s
   367   331   62.99272   12   55   50.00000   62.99272  26.0%  47.4  432s
   375   336   62.96010   13   44   50.00000   62.99272  26.0%  50.8  436s
   383   341   62.86096   14   42   50.00000   62.99272  26.0%  53.3  440s
   395   349   62.11111   15   15   50.00000   62.99272  26.0%  54.9  447s
   403   356   61.25000   16   34   50.00000   62.99272  26.0%  56.1  450s
   412   361   62.40000   17   20   50.00000   62.99272  26.0%  59.9  455s
   416   363   62.40000   18    8   50.00000   62.99272  26.0%  61.1  461s
   420   366   61.31286   18   43   50.00000   62.99272  26.0%  61.9  476s
   429   374   62.35294   19   13   50.00000   62.99272  26.0%  62.9  480s
   442   384   62.29723   21   50   50.00000   62.99272  26.0%  63.8  485s
   458   396   62.29275   23   50   50.00000   62.99272  26.0%  69.5  492s
   467   403   62.28421   24   16   50.00000   62.99272  26.0%  71.6  495s
   489   411   59.28889   26   24   50.00000   62.99272  26.0%  81.3  503s
   497   424   62.27544   27   33   50.00000   62.99272  26.0%  83.8  508s
   513   430   62.26671   29   29   50.00000   62.99272  26.0%  87.0  512s
   527   435   62.26296   30   23   50.00000   62.99272  26.0%  91.9  515s
   562   452   62.25516   33   45   50.00000   62.99272  26.0%  94.7  523s
   577   464   62.25351   35   32   50.00000   62.99272  26.0%   102  529s
   594   472   62.25029   37   27   50.00000   62.99272  26.0%   107  534s
   608   481   59.98667   38   39   50.00000   62.99272  26.0%   112  539s
   624   488   62.24442   41   35   50.00000   62.99272  26.0%   115  543s
   638   493   60.66667   43   13   50.00000   62.99272  26.0%   124  548s
   658   505   62.23658   44   46   50.00000   62.99272  26.0%   123  554s
   686   517   62.22725   46   45   50.00000   62.99272  26.0%   123  559s
   708   528   62.22636   47   36   50.00000   62.99272  26.0%   126  564s
   726   536   62.20935   50   56   50.00000   62.99272  26.0%   130  569s
   746   546   62.20697   52   56   50.00000   62.99272  26.0%   135  576s
   765   560   62.19752   55   29   50.00000   62.99272  26.0%   139  584s
   789   570   62.19413   57   35   50.00000   62.99272  26.0%   143  588s
   813   565     cutoff   59        50.00000   62.99272  26.0%   146  604s
H  814   534                      54.0000000   62.99272  16.7%   146  604s
H  815   363                      59.0000000   62.99272  6.77%   146  604s
H  816   348                      60.0000000   62.99272  4.99%   146  604s
   818   344     cutoff   60        60.00000   62.99272  4.99%   146  608s
   882   338   62.18854   64   34   60.00000   62.99272  4.99%   145  613s
   951   326   62.95358   24   42   60.00000   62.99272  4.99%   143  618s
* 1009   293              27      61.0000000   62.98422  3.25%   139  618s
  1055   240   62.67102   19   39   61.00000   62.92606  3.16%   135  624s
  1132   197 infeasible   18        61.00000   62.81818  2.98%   134  632s
  1230   177 infeasible   21        61.00000   62.73681  2.85%   130  639s
  1345    51   62.31475   21   35   61.00000   62.60114  2.62%   123  645s

Cutting planes:
  Gomory: 3
  Cover: 1
  Implied bound: 19
  MIR: 21
  Mixing: 1
  StrongCG: 1
  Flow cover: 94
  Inf proof: 2
  Zero half: 1
  Mod-K: 1
  RLT: 1

Explored 1577 nodes (175523 simplex iterations) in 647.07 seconds (337.60 work units)
Thread count was 8 (of 8 available processors)

Solution count 10: 61 60 59 ... 7

Optimal solution found (tolerance 1.00e-04)
Best objective 6.100000000000e+01, best bound 6.100000000000e+01, gap 0.0000%
最佳路径: [(0, 211), (4, 100), (21, 150), (23, 273), (25, 171), (34, 126), (44, 280), (48, 162), (49, 215), (57, 208), (64, 170), (65, 114), (69, 173), (73, 132), (75, 73), (78, 116), (88, 21), (92, 276), (96, 0), (100, 175), (105, 161), (106, 25), (114, 78), (116, 154), (119, 260), (122, 23), (126, 57), (132, 298), (134, 135), (135, 287), (143, 232), (150, 241), (154, 259), (158, 64), (161, 119), (162, 198), (170, 34), (171, 192), (173, 214), (175, 184), (184, 222), (192, 290), (198, 256), (208, 288), (211, 49), (214, 158), (215, 48), (222, 96), (232, 69), (241, 92), (256, 122), (259, 75), (260, 44), (266, 106), (273, 65), (276, 266), (280, 143), (287, 105), (288, 4), (290, 134), (298, 88)]
花费时间: 142.5
访问的城市数量: 61

考虑到300阶过大，换成50阶的，用时0.1s：

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
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (win64 - Windows 11+.0 (26120.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 2551 rows, 2550 columns and 19259 nonzeros
Model fingerprint: 0x1d9ad627
Variable types: 50 continuous, 2500 integer (2500 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+03]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+02]
Presolve removed 1432 rows and 2438 columns
Presolve time: 0.08s
Presolved: 1119 rows, 112 columns, 2786 nonzeros
Variable types: 33 continuous, 79 integer (79 binary)

Root relaxation: objective 1.276923e+01, 175 iterations, 0.00 seconds (0.00 work units)

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

     0     0   12.76923    0   10          -   12.76923      -     -    0s
     0     0   12.47012    0   10          -   12.47012      -     -    0s
H    0     0                       5.0000000   12.47012   149%     -    0s
     0     0   11.37942    0   16    5.00000   11.37942   128%     -    0s
H    0     0                       7.0000000   11.25101  60.7%     -    0s
H    0     0                       9.0000000   11.02273  22.5%     -    0s
     0     0   11.02273    0   16    9.00000   11.02273  22.5%     -    0s
     0     0   11.02273    0   10    9.00000   11.02273  22.5%     -    0s
H    0     0                      10.0000000   11.02273  10.2%     -    0s

Cutting planes:
  Gomory: 1
  Implied bound: 6

Explored 1 nodes (210 simplex iterations) in 0.20 seconds (0.10 work units)
Thread count was 8 (of 8 available processors)

Solution count 4: 10 9 7 5 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0000%
最佳路径: [(0, 23), (1, 12), (4, 46), (12, 4), (13, 35), (19, 1), (23, 13), (35, 19), (42, 0), (46, 42)]
花费时间: 120.5
访问的城市数量: 10

3.2 Matlab 求解器

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

% 用时间矩阵求解最大城市数目
disp('求解开始时间为：');
s = datetime('now')
Gmatrix = randi([1, 1000], 300, 300);
max_time = 144 - 2.5; % 144减去在广州的2.5小时
[path, total_time, cities_visited] = solve_max_cities_with_time_constraint(Gmatrix, max_time, 4);
disp('最佳路径:');
disp(path);
disp('花费时间:');
disp(total_time + 2.5);
disp('访问的城市数量:');
disp(cities_visited);
disp('求解结束时间为：');
e = datetime('now')
disp('求解耗时：' , e - s);


function [path, total_time, cities_visited] = solve_max_cities_with_time_constraint(Gmatrix, max_time, start_city)
    n = size(Gmatrix, 1);
    % 输入参数：时间矩阵 Gmatrix，最大时间 max_time，起点城市 start_city
    % 输出参数：最佳路径 path，花费时间 total_time，访问的城市数量 cities_visited

    % 创建模型
    model = optimproblem('ObjectiveSense', 'maximize');

    % 决策变量：x(i, j) ,是否选择从城市 i 到城市 j 的路径
    x = optimvar('x', n, n, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);

    % 辅助变量：u(i) 用于避免子环
    u = optimvar('u', n, 'Type', 'continuous', 'LowerBound', 0);

    % 最大化访问的城市数量
    model.Objective = sum(x, 'all');

    % 总时间花费 <= max_time
    model.Constraints.time_limit = sum(sum(Gmatrix .* x)) <= max_time;

    % 添加流量约束
    for i = 1:n-1
        model.Constraints.(['flow_balance_', num2str(i)]) = sum(x(i, :)) == sum(x(:, i));
        model.Constraints.(['out_', num2str(i)]) = sum(x(i, :)) <= 1;
        model.Constraints.(['in_', num2str(i)]) = sum(x(:, i)) <= 1;
    end

    % 添加起点约束：确保从指定的起点城市出发
    model.Constraints.start_city_out = sum(x(start_city, :)) == 1;

    % 添加子环消除约束（MTZ约束）
    for i = 2:n
        for j = 2:n
            if i ~= j
                model.Constraints.(['mtz_', num2str(i), '_', num2str(j)]) = u(i) - u(j) + n * x(i, j) <= n - 1;
            end
        end
    end

    % 防止自环
    for i = 1:n
        model.Constraints.(['no_self_loop_', num2str(i)]) = x(i, i) == 0;
    end

    % 求解模型
    options = optimoptions('intlinprog', 'Display', 'iter');
    [sol, fval, exitflag] = solve(model, 'Options', options);

    % 检查是否找到最优解
    if exitflag == 1
        path = [];
        for i = 1:n
            for j = 1:n
                if sol.x(i, j) > 0.5  % 判断变量值是否为1
                    path = [path; i, j];
                end
            end
        end
        total_time = sum(sum(Gmatrix .* sol.x));
        cities_visited = length(unique([path(:, 1); path(:, 2)]));
    else
        disp('没有找到最优解');
        path = [];
        total_time = 0;
        cities_visited = 0;
    end
end

遗憾的是在300阶的情况下，Matlab求解了1个小时都没有求解出来，可能永远也求不出来了。

换成50阶的,用时124s：

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
>> solve
求解开始时间为：

s = 

  datetime

   2024-08-31 18:23:42



将使用 intlinprog 求解问题。

LP:                Optimal objective value is 10.651563.                                            


Cut Generation:    Applied 3 Gomory cuts, 5 mir cuts,                                               
                   and 1 cover cut.                                                                 
                   Upper bound is 9.000000.                                                         


Branch and Bound:


   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
      24      0.68         1   7.000000e+00   1.250000e+01                                          
      27      0.69         1   7.000000e+00   1.250000e+01                                          

找到最优解。

Intlinprog 已停止，因为目标值在最佳值的间隙容差范围内，options.RelativeGapTolerance = 0.0001。intcon 变量是容差范围内的整
数，options.IntegerTolerance = 1e-05。


最佳路径:
     1    44
     4    13
    13     1
    27    37
    29    27
    37     4
    44    29


花费时间:
  124.5000


访问的城市数量:
     7


求解结束时间为：

e = 

  datetime

   2024-08-31 18:24:23

#math,coding
TSP问题的优化求解
https://silenzio111.github.io/2024/08/31/opttsp/
作者
silenzio
发布于
2024年8月31日
许可协议