from gurobipy import Model, quicksum, GRB
#from mypulp import Model, quicksum, GRB
#from mindoptpy import Model, quicksum, Column
#from mindoptpy import MDO as GRB
import math
import networkx as nx時間枠付き巡回セールスマン問題
時間枠付き巡回セールスマン問題に対する定式化とアルゴリズム
準備
時間枠付き巡回セールスマン問題
ここでは，巡回セールスマン問題に時間枠を追加した 時間枠付き巡回セールスマン問題(traveling salesman problem with time windows)を考える．
この問題は，特定の点 \(0\) を時刻 \(0\) に出発すると仮定し， 点間の移動距離 \(c_{ij}\) を移動時間とみなし， さらに点 \(i\) に対する出発時刻が最早時刻 \(e_i\) と最遅時刻 \(\ell_i\) の間でなければならないという制約を課した問題である． ただし，時刻 \(e_i\) より早く点 \(i\) に到着した場合には，点 \(i\) 上で時刻 \(e_i\) まで待つことができるものとする．
ポテンシャル定式化
巡回セールスマン問題に対するポテンシャル制約の拡張を考える．
点 \(i\) を出発する時刻を表す変数 \(t_i\) を導入する． \(t_i\) は以下の制約を満たす必要がある． \[ e_i \leq t_i \leq \ell_i \ \ \ \forall i=1,2,\ldots,n \] ただし， \(e_1=0, \ell_1=\infty\) と仮定する．
点 \(i\) の次に点 \(j\) を訪問する \((x_{ij}=1)\) ときには， 点 \(j\) を出発する時刻 \(t_j\) は，点 \(i\) を出発する時刻に移動時間 \(c_{ij}\) を加えた値以上であることから， 以下の式を得る． \[ t_i + c_{ij} - M (1-x_{ij}) \leq t_j \ \ \ \forall i,j, j \neq 1, i \neq j \] ここで，\(M\) は大きな数を表す定数である． なお，移動時間 \(c_{ij}\) は正の数と仮定する．\(c_{ij}\) が \(0\) だと \(t_i=t_j\) になる可能性があり， 部分巡回路ができてしまう．これを避けるためには，巡回セールスマン問題と同様の制約を付加する必要があるが， \(c_{ij}>0\) の仮定の下では，上の制約によって部分巡回路を除去することができる．
このような大きな数Big Mを含んだ定式化はあまり実用的ではないので，時間枠を用いて強化したものを示す．
\[ \begin{array}{lll} minimize & \sum_{i \neq j} c_{ij} x_{ij} & \\ s.t. & \sum_{j: j \neq i} x_{ij} = 1 & \forall i=1,2,\ldots,n \\ & \sum_{j: j \neq i} x_{ji} = 1 & \forall i=1,2,\ldots,n \\ & t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) \leq t_j & \forall i,j,j \neq 1, i \neq j \\ & x_{ij} \in \{0,1\} & \forall i,j, i \neq j \\ & e_i \leq t_{i} \leq \ell_i & \forall i=1,2,\ldots,n \end{array} \]
巡回セールスマン問題のときと同様に，ポテンシャル制約と上下限制約は， 持ち上げ操作によってさらに以下のように強化できる．
\[ t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) + [\ell_i-e_j+\min \{-c_{ji}, e_j-e_i\} ]^+ x_{ji} \leq t_j \ \ \ \forall i,j, j \neq1, i \neq j \]
\[ e_i + \sum_{j \neq i} [e_j+ c_{ji} -e_i]^{+} x_{ji} \leq t_i \ \ \ \forall i=2,\cdots,n \]
\[ t_i \leq \ell_i -\sum_{j \neq 1,i} [\ell_i -\ell_j+c_{ij}]^+ x_{ij} \ \ \ \forall i=2,\cdots,n \]
以下に，強化していない標準定式化 mtztw と強化した定式化 mtz2tw のコードを示す．
def mtztw(n, c, e, l):
"""mtzts: model for the traveling salesman problem with time windows
(based on Miller-Tucker-Zemlin's one-index potential formulation)
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
- e[i]: earliest date for visiting node i
- l[i]: latest date for visiting node i
Returns a model, ready to be solved.
"""
model = Model("tsptw - mtz")
x, u, v = {}, {}, {}
for i in range(1, n + 1):
u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
v[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name="v(%s)" % i)
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n + 1):
model.addConstr(
quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
)
for i in range(1, n + 1):
for j in range(2, n + 1):
if i != j:
M = max(l[i] + c[i, j] - e[j], 0)
model.addConstr(
u[i] - u[j] + M * x[i, j] <= M - c[i, j], "MTZ(%s,%s)" % (i, j)
)
model.addConstr(
v[i] - v[j] + (n - 1) * x[i, j] <= n - 2, "MTZorig(%s,%s)" % (i, j)
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.update()
model.__data = x, u
return model
def mtz2tw(n, c, e, l):
"""mtz: model for the traveling salesman problem with time windows
(based on Miller-Tucker-Zemlin's one-index potential formulation, stronger constraints)
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
- e[i]: earliest date for visiting node i
- l[i]: latest date for visiting node i
Returns a model, ready to be solved.
"""
model = Model("tsptw - mtz-strong")
x, u = {}, {}
for i in range(1, n + 1):
u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n + 1):
model.addConstr(
quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
)
for j in range(2, n + 1):
if i != j:
M1 = max(l[i] + c[i, j] - e[j], 0)
M2 = max(l[i] + min(-c[j, i], e[j] - e[i]) - e[j], 0)
model.addConstr(
u[i] + c[i, j] - M1 * (1 - x[i, j]) + M2 * x[j, i] <= u[j],
"LiftedMTZ(%s,%s)" % (i, j),
)
for i in range(2, n + 1):
model.addConstr(
e[i]
+ quicksum(
max(e[j] + c[j, i] - e[i], 0) * x[j, i]
for j in range(1, n + 1)
if i != j
)
<= u[i],
"LiftedLB(%s)" % i,
)
model.addConstr(
u[i]
<= l[i]
- quicksum(
max(l[i] - l[j] + c[i, j], 0) * x[i, j]
for j in range(2, n + 1)
if i != j
),
"LiftedUB(%s)" % i,
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.update()
model.__data = x, u
return model# hide
def distance(x1, y1, x2, y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
def make_data(n, width):
"""make_data: compute matrix distance and time windows."""
x = dict([(i, 100 * random.random()) for i in range(1, n + 1)])
y = dict([(i, 100 * random.random()) for i in range(1, n + 1)])
c = {}
for i in range(1, n + 1):
for j in range(1, n + 1):
if j != i:
c[i, j] = distance(x[i], y[i], x[j], y[j])
e = {1: 0}
l = {1: 0}
start = 0
delta = int(76.0 * math.sqrt(n) / n * width) + 1
for i in range(1, n):
j = i + 1
start += c[i, j]
e[j] = max(start - delta, 0)
l[j] = start + delta
return c, x, y, e, l
# this formulation is too slow!
def tsptw2(n, c, e, l):
"""tsptw2: model for the traveling salesman problem with time windows
(based on Miller-Tucker-Zemlin's formulation, two-index potential)
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
- e[i]: earliest date for visiting node i
- l[i]: latest date for visiting node i
Returns a model, ready to be solved.
"""
model = Model("tsptw2")
x, u = {}, {}
for i in range(1, n + 1):
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
u[i, j] = model.addVar(vtype="C", name="u(%s,%s)" % (i, j))
model.update()
for i in range(1, n + 1):
model.addConstr(
quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
)
for j in range(2, n + 1):
model.addConstr(
quicksum(u[i, j] + c[i, j] * x[i, j] for i in range(1, n + 1) if i != j)
- quicksum(u[j, k] for k in range(1, n + 1) if k != j)
<= 0,
"Relate(%s)" % j,
)
for i in range(1, n + 1):
for j in range(1, n + 1):
if i != j:
model.addConstr(e[i] * x[i, j] <= u[i, j], "LB(%s,%s)" % (i, j))
model.addConstr(u[i, j] <= l[i] * x[i, j], "UB(%s,%s)" % (i, j))
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.update()
model.__data = x, u
return modelベンチマーク問題例の読み込み
時間枠付き巡回セールスマン問題のベンチマーク問題例は，以下のサイトから入手できる．
- 200点までの古典的なベンチマーク問題例: https://lopez-ibanez.eu/tsptw-instances
- 400点までの新しいベンチマーク問題例： https://homepages.dcc.ufmg.br/~rfsilva/tsptw/#tools
読み込んでデータを準備する．
import math
folder = "../data/tsptw/"
fn = "n60w20.001.txt"
#fn ="n200w20.002.txt" #116s
# fn ="n100w40.002.txt"
# fn ="n100w60.002.txt"
#fn ="n250w400.001.txt" #191s
#fn = "n350w100.001.txt" #5.84s
#fn ="n400w500.005.txt" #FEASIBLE 12843
#fn = "n400w100.001.txt" #19s
#fn = "n400w300.001.txt" #70s
#fn = "n400w400.001.txt"
#fn ="n200w40.001.txt"
f = open(folder + fn)
data = f.readlines()
f.close()x_, y_, early, late = [], [], [], []
for row in data:
try:
L = list(map(float, row.split()))
if L[0] >= 900:
continue
x_.append(L[1])
y_.append(L[2])
early.append(L[4])
late.append(L[5])
except:
pass
n = len(x_)
e, l = {}, {}
for i in range(n):
e[i + 1] = int(early[i])
l[i + 1] = int(late[i])
c = {}
for i in range(n):
for j in range(n):
c[i + 1, j + 1] = int(math.sqrt((x_[i] - x_[j]) ** 2 + (y_[i] - y_[j]) ** 2))for i in range(1, n + 1):
for j in range(1, n + 1):
for k in range(1, n + 1):
c[i, j] = min(c[i, j], c[i, k] + c[k, j])求解と描画
まず，強化していない標準定式化を用いて求解する．
model = mtztw(n, c, e, l)
model.optimize()
x, u = model.__data
sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]
print("Opt.value =", model.ObjVal, sol)次に，時間枠を用いて制約を強化した定式化を用いて求解する．
model = mtz2tw(n, c, e, l)
model.optimize()
x, u = model.__data
sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]
print("Opt.value =", model.ObjVal, sol)Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (mac64[x86])
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3842 rows, 3721 columns and 22227 nonzeros
Model fingerprint: 0x10540bf3
Variable types: 61 continuous, 3660 integer (3660 binary)
Coefficient statistics:
Matrix range [1e+00, 6e+02]
Objective range [1e+00, 6e+01]
Bounds range [1e+00, 6e+02]
RHS range [1e+00, 6e+02]
Presolve removed 2118 rows and 2215 columns
Presolve time: 0.05s
Presolved: 1724 rows, 1506 columns, 7938 nonzeros
Variable types: 57 continuous, 1449 integer (1449 binary)
Root relaxation: objective 5.361061e+02, 336 iterations, 0.01 seconds (0.01 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 536.10606 0 67 - 536.10606 - - 0s
0 0 541.00000 0 13 - 541.00000 - - 0s
H 0 0 553.0000000 541.00000 2.17% - 0s
0 0 541.00000 0 10 553.00000 541.00000 2.17% - 0s
H 0 0 552.0000000 541.00000 1.99% - 0s
H 0 0 551.0000000 541.00000 1.81% - 0s
0 0 541.00000 0 6 551.00000 541.00000 1.81% - 0s
0 0 541.00000 0 42 551.00000 541.00000 1.81% - 0s
0 0 541.00000 0 16 551.00000 541.00000 1.81% - 0s
0 0 541.05935 0 11 551.00000 541.05935 1.80% - 0s
0 0 541.05935 0 11 551.00000 541.05935 1.80% - 0s
0 0 541.37202 0 11 551.00000 541.37202 1.75% - 0s
0 0 541.42378 0 11 551.00000 541.42378 1.74% - 0s
0 0 541.76983 0 16 551.00000 541.76983 1.68% - 0s
0 0 541.78378 0 18 551.00000 541.78378 1.67% - 0s
0 0 542.14093 0 18 551.00000 542.14093 1.61% - 0s
0 0 542.21128 0 17 551.00000 542.21128 1.60% - 0s
0 0 542.22222 0 26 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 11 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 13 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 11 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 11 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 11 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 11 551.00000 542.22222 1.59% - 0s
0 0 542.22222 0 12 551.00000 542.22222 1.59% - 0s
0 0 542.50824 0 11 551.00000 542.50824 1.54% - 0s
0 0 542.50824 0 10 551.00000 542.50824 1.54% - 0s
0 2 542.50824 0 10 551.00000 542.50824 1.54% - 0s
Cutting planes:
Learned: 5
Gomory: 5
Cover: 2
MIR: 4
StrongCG: 2
GUB cover: 2
Zero half: 5
RLT: 3
Relax-and-lift: 3
Explored 348 nodes (3140 simplex iterations) in 0.48 seconds (0.26 work units)
Thread count was 16 (of 16 available processors)
Solution count 3: 551 552 553
Optimal solution found (tolerance 1.00e-04)
Best objective 5.510000000000e+02, best bound 5.510000000000e+02, gap 0.0000%
Opt.value = 551.0 [1, 7, 13, 39, 49, 10, 16, 52, 35, 57, 18, 20, 43, 59, 23, 33, 60, 8, 26, 14, 5, 53, 6, 45, 51, 31, 48, 19, 47, 32, 17, 24, 22, 3, 29, 42, 11, 30, 12, 28, 40, 44, 46, 36, 27, 9, 58, 4, 34, 55, 2, 25, 50, 38, 21, 37, 56, 15, 41, 61, 54]G = nx.Graph()
for idx, i in enumerate(sol[:-1]):
G.add_edge(i, sol[idx + 1])
G.add_edge(sol[-1],sol[1])
pos = {i: (x_[i-1], y_[i-1]) for i in range(1,n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
OR-Tools(cp-sat)による求解
ここでは，OR-Tools(cp-sat)を用いて時間枠付き巡回セールスマン問題を解く．
巡回路に含まれる枝 \((i,j)\) を論理変数 \(x_{ij}\) で表し，各点 \(i\) への到着時刻を整数変数 \(u_i\) で表す．
cp-satでは， 変数 \(x\) が巡回路になっていることを表現するための閉路制約が準備されている． モデルインスタンスのAddCircuitメソッドで制約を定義する．引数は， 枝の始点 \(i\) ，終点 \(j\)， 論理変数 \(x_{ij}\) ３つ組（タプル）のリストである． 巡回路に含まれない点 \(i\) に対しては， 自己閉路を表す \(x_{ii}\) がTrueになる．
到着時刻を表す変数 \(u\) と枝変数 \(x\) を繋ぐためには，論理制約を使うことができる． これは，モデルに線形制約をAddメソッドで追加する際に， 追加した時間差を表す制約式の後にOnlyEnforceIfメソッドを付加することによって定義される． OnlyEnforceIfメソッドの引数は論理変数であり，この場合には\(x_{ij}\)となる． つまり， \(x_{ij}\) がTrueのときに \(u_i+c_{ij} \leq u_j\) が付加されることになる．
from ortools.sat.python import cp_model
model = cp_model.CpModel()
x, u = {}, {}
for i in range(1, n + 1):
u[i] = model.NewIntVar(lb=e[i], ub=l[i], name=f"u({i})")
for j in range(1, n + 1):
if i != j:
x[i, j] = model.NewBoolVar(name=f"x({i},{j})")
circuit = [(i,j,x_) for (i, j), x_ in x.items()]
model.AddCircuit(circuit) #枝集合が閉路になっている制約を表す
model.Add(u[1] == 0) # fix the root node to depth 0
for i in range(1, n + 1):
for j in range(1, n + 1):
if i != j:
if j != 1: # The root node is special.
model.Add(u[i] + c[i,j] <= u[j]).OnlyEnforceIf(x[i, j])
model.Minimize(sum(c[i,j] * x_ for (i,j), x_ in x.items()))
solver = cp_model.CpSolver()
#solver.parameters.max_time_in_seconds = 10
solver.parameters.log_search_progress = True
status = solver.Solve(model)
print("status", solver.StatusName(status))
Starting CP-SAT solver v9.9.3963
Parameters: log_search_progress: true
Setting number of workers to 12
Initial optimization model '': (model_fingerprint: 0x4d053609fffe6e2f)
#Variables: 3'721 (#bools: 3'660 in objective)
- 3'660 Booleans in [0,1]
- 1 in [0,23]
- 1 in [0,27]
- 1 in [0,601]
- 1 in [4,31]
- 1 in [7,43]
- 1 in [11,18]
- 1 in [16,36]
- 1 in [23,48]
- 1 in [24,38]
- 1 in [26,47]
- 1 in [30,51]
- 1 in [32,41]
- 1 in [52,77]
- 1 in [67,70]
- 1 in [68,101]
- 1 in [71,90]
- 1 in [73,85]
- 1 in [79,104]
- 1 in [105,131]
- 1 in [107,112]
- 1 in [118,141]
- 1 in [140,146]
- 1 in [148,160]
- 1 in [156,180]
- 1 in [164,177]
- 1 in [180,199]
- 1 in [185,203]
- 1 in [193,215]
- 1 in [195,213]
- 1 in [200,214]
- 1 in [207,237]
- 1 in [211,227]
- 1 in [222,248]
- 1 in [227,248]
- 1 in [234,245]
- 1 in [237,260]
- 1 in [238,257]
- 1 in [251,278]
- 1 in [258,261]
- 1 in [263,284]
- 1 in [273,300]
- 1 in [291,308]
- 1 in [309,322]
- 1 in [314,340]
- 1 in [318,334]
- 1 in [327,338]
- 1 in [327,360]
- 1 in [333,365]
- 1 in [350,361]
- 1 in [360,368]
- 1 in [371,376]
- 1 in [387,405]
- 1 in [395,410]
- 1 in [401,421]
- 1 in [405,420]
- 1 in [435,444]
- 1 in [452,481]
- 1 in [473,509]
- 1 in [525,541]
- 1 in [552,583]
- 1 in [559,583]
#kCircuit: 1
#kLinear1: 1
#kLinear2: 3'600 (#enforced: 3'600)
Starting presolve at 0.00s
9.33e-04s 0.00e+00d [DetectDominanceRelations]
3.43e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=6 #num_dual_strengthening=5
8.00e-06s 0.00e+00d [ExtractEncodingFromLinear]
1.11e-02s 1.73e-02d [Probe] #probed=3'549 #fixed_bools=365 #new_bounds=5 #equiv=3 #new_binary_clauses=321
1.10e-05s 0.00e+00d [MaxClique]
7.02e-04s 0.00e+00d [DetectDominanceRelations]
1.67e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=3 #num_dual_strengthening=1
3.90e-05s 0.00e+00d [ProcessAtMostOneAndLinear]
3.53e-04s 0.00e+00d [DetectDuplicateConstraints]
1.50e-05s 1.80e-08d [DetectDominatedLinearConstraints] #relevant_constraints=3
2.10e-04s 0.00e+00d [DetectDifferentVariables] #different=4
1.00e-05s 0.00e+00d [ProcessSetPPC]
9.00e-06s 0.00e+00d [FindAlmostIdenticalLinearConstraints]
6.20e-05s 2.60e-06d [FindBigAtMostOneAndLinearOverlap]
6.00e-05s 3.60e-05d [FindBigVerticalLinearOverlap]
7.00e-06s 0.00e+00d [FindBigHorizontalLinearOverlap]
3.20e-04s 0.00e+00d [MergeClauses]
4.38e-04s 0.00e+00d [DetectDominanceRelations]
8.74e-04s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
6.82e-04s 0.00e+00d [DetectDominanceRelations]
1.15e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
[SAT presolve] num removable Booleans: 0 / 1595
[SAT presolve] num trivial clauses: 0
[SAT presolve] [0s] clauses:65 literals:130 vars:130 one_side_vars:130 simple_definition:0 singleton_clauses:0
[SAT presolve] [9e-06s] clauses:65 literals:130 vars:130 one_side_vars:130 simple_definition:0 singleton_clauses:0
[SAT presolve] [4.8e-05s] clauses:65 literals:130 vars:130 one_side_vars:130 simple_definition:0 singleton_clauses:0
9.12e-03s 1.42e-02d [Probe] #probed=3'055 #fixed_bools=189 #new_bounds=3 #equiv=3 #new_binary_clauses=196
1.51e-04s 0.00e+00d [MaxClique] Merged 65(130 literals) into 63(242 literals) at_most_ones.
7.03e-04s 0.00e+00d [DetectDominanceRelations]
1.45e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=3 #num_dual_strengthening=1
7.00e-05s 0.00e+00d [ProcessAtMostOneAndLinear]
3.44e-04s 0.00e+00d [DetectDuplicateConstraints]
3.25e-04s 2.40e-08d [DetectDominatedLinearConstraints] #relevant_constraints=4
2.10e-04s 0.00e+00d [DetectDifferentVariables] #different=7
4.10e-05s 7.16e-07d [ProcessSetPPC] #relevant_constraints=59
3.22e-04s 0.00e+00d [FindAlmostIdenticalLinearConstraints]
8.00e-05s 2.30e-05d [FindBigAtMostOneAndLinearOverlap]
6.00e-05s 3.49e-05d [FindBigVerticalLinearOverlap]
7.00e-06s 0.00e+00d [FindBigHorizontalLinearOverlap]
9.00e-06s 0.00e+00d [MergeClauses]
6.87e-04s 0.00e+00d [DetectDominanceRelations]
1.11e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
3.71e-04s 0.00e+00d [DetectDominanceRelations]
7.80e-04s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
[SAT presolve] num removable Booleans: 0 / 1404
[SAT presolve] num trivial clauses: 0
[SAT presolve] [0s] clauses:61 literals:122 vars:122 one_side_vars:122 simple_definition:0 singleton_clauses:0
[SAT presolve] [8e-06s] clauses:61 literals:122 vars:122 one_side_vars:122 simple_definition:0 singleton_clauses:0
[SAT presolve] [4.6e-05s] clauses:61 literals:122 vars:122 one_side_vars:122 simple_definition:0 singleton_clauses:0
7.63e-03s 1.10e-02d [Probe] #probed=2'607 #fixed_bools=151 #equiv=4 #new_binary_clauses=118
1.57e-04s 0.00e+00d [MaxClique] Merged 112(340 literals) into 59(234 literals) at_most_ones.
3.91e-04s 0.00e+00d [DetectDominanceRelations]
8.97e-04s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
6.20e-05s 0.00e+00d [ProcessAtMostOneAndLinear]
3.20e-04s 0.00e+00d [DetectDuplicateConstraints]
3.29e-04s 2.40e-08d [DetectDominatedLinearConstraints] #relevant_constraints=4
2.05e-04s 0.00e+00d [DetectDifferentVariables] #different=7
3.80e-05s 7.16e-07d [ProcessSetPPC] #relevant_constraints=59
3.28e-04s 0.00e+00d [FindAlmostIdenticalLinearConstraints]
7.90e-05s 2.29e-05d [FindBigAtMostOneAndLinearOverlap]
6.20e-05s 3.32e-05d [FindBigVerticalLinearOverlap]
7.00e-06s 0.00e+00d [FindBigHorizontalLinearOverlap]
8.00e-06s 0.00e+00d [MergeClauses]
7.09e-04s 0.00e+00d [DetectDominanceRelations]
1.14e-03s 0.00e+00d [PresolveToFixPoint] #num_loops=1 #num_dual_strengthening=1
1.79e-04s 0.00e+00d [ExpandObjective]
Presolve summary:
- 13 affine relations were detected.
- rule 'TODO dual: add implied bound' was applied 146 times.
- rule 'TODO dual: only one blocking constraint?' was applied 21'296 times.
- rule 'TODO dual: only one blocking enforced constraint?' was applied 146 times.
- rule 'TODO dual: only one unspecified blocking constraint?' was applied 6 times.
- rule 'affine: new relation' was applied 13 times.
- rule 'at_most_one: removed literals' was applied 2 times.
- rule 'at_most_one: size one' was applied 2 times.
- rule 'at_most_one: transformed into max clique.' was applied 2 times.
- rule 'bool_or: only one literal' was applied 1'697 times.
- rule 'circuit: degree 2' was applied 3 times.
- rule 'circuit: removed false arcs.' was applied 4 times.
- rule 'deductions: 268 stored' was applied 1 time.
- rule 'dual: reduced domain' was applied 27 times.
- rule 'enforcement: false literal' was applied 11 times.
- rule 'enforcement: true literal' was applied 4 times.
- rule 'incompatible linear: add implication' was applied 194 times.
- rule 'linear2: implied ax + by = cte has only one solution' was applied 5 times.
- rule 'linear: always true' was applied 1'629 times.
- rule 'linear: empty' was applied 1 time.
- rule 'linear: fixed or dup variables' was applied 61 times.
- rule 'linear: reduced variable domains' was applied 1 time.
- rule 'linear: simplified rhs' was applied 387 times.
- rule 'objective: variable not used elsewhere' was applied 2'398 times.
- rule 'presolve: 2399 unused variables removed.' was applied 1 time.
- rule 'presolve: iteration' was applied 3 times.
- rule 'variables: detect half reified value encoding' was applied 10 times.
Presolved optimization model '': (model_fingerprint: 0xd9b43925f5d11727)
#Variables: 1'306 (#bools: 1'245 in objective)
- 1'245 Booleans in [0,1]
- 1 in [4,31]
- 1 in [5,23]
- 1 in [6,27]
- 1 in [11,18]
- 1 in [14,43]
- 1 in [16,36]
- 1 in [23,48]
- 1 in [24,38]
- 1 in [26,47]
- 1 in [30,51]
- 1 in [32,41]
- 1 in [55,61]
- 1 in [67,70]
- 1 in [68,101]
- 1 in [71,90]
- 1 in [73,85]
- 1 in [79,104]
- 1 in [109,112]
- 1 in [118,121]
- 1 in [135,138]
- 1 in [141,144]
- 1 in [148,151]
- 1 in [157,180]
- 1 in [164,177]
- 1 in [180,188]
- 1 in [189,203]
- 1 in [193,215]
- 1 in [195,213]
- 1 in [200,214]
- 1 in [207,237]
- 1 in [211,227]
- 1 in [222,248]
- 1 in [227,248]
- 1 in [234,245]
- 1 in [237,260]
- 1 in [238,257]
- 1 in [251,278]
- 1 in [258,261]
- 1 in [263,284]
- 1 in [284,296]
- 1 in [297,308]
- 1 in [309,322]
- 1 in [314,340]
- 1 in [318,334]
- 1 in [327,338]
- 1 in [327,360]
- 1 in [333,365]
- 1 in [350,361]
- 1 in [362,367]
- 1 in [371,376]
- 1 in [389,405]
- 1 in [395,410]
- 1 in [402,409][415,421]
- 1 in [405,419]
- 1 in [441,444]
- 1 in [467,470]
- 1 in [488,491]
- 1 in [525,528]
- 1 in [561,583]
- 1 in [568,583]
- 1 constants in {1}
#kAtMostOne: 51 (#literals: 218)
#kBoolAnd: 69 (#enforced: 69) (#literals: 138)
#kCircuit: 1
#kLinear1: 8 (#enforced: 8)
#kLinear2: 260 (#enforced: 256)
Preloading model.
#Bound 0.05s best:inf next:[212,29336] initial_domain
#Model 0.05s var:1305/1306 constraints:389/389
Starting search at 0.05s with 12 workers.
8 full problem subsolvers: [core, default_lp, max_lp, no_lp, pseudo_costs, quick_restart, quick_restart_no_lp, reduced_costs]
3 first solution subsolvers: [fj_long_default, fj_short_default, fs_random]
11 incomplete subsolvers: [feasibility_pump, graph_arc_lns, graph_cst_lns, graph_dec_lns, graph_var_lns, rins/rens, rnd_cst_lns, rnd_var_lns, routing_path_lns, routing_random_lns, violation_ls]
3 helper subsolvers: [neighborhood_helper, synchronization_agent, update_gap_integral]
#Bound 0.06s best:inf next:[212,27763] core (initial_propagation)
#Bound 0.06s best:inf next:[520,27763] quick_restart (initial_propagation)
#Bound 0.06s best:inf next:[541,27763] max_lp
#Model 0.07s var:1217/1306 constraints:383/389
#1 0.07s best:552 next:[541,551] default_lp (fixed_bools=90/1273)
#Bound 0.07s best:552 next:[542,551] max_lp
#Model 0.07s var:630/1306 constraints:380/389
#Bound 0.08s best:552 next:[543,551] max_lp
#Model 0.09s var:628/1306 constraints:380/389
#Model 0.09s var:619/1306 constraints:375/389
#Model 0.10s var:506/1306 constraints:358/389
#Bound 0.11s best:552 next:[544,551] pseudo_costs
#2 0.11s best:551 next:[544,550] default_lp (fixed_bools=1009/1278)
#Model 0.11s var:274/1306 constraints:311/389
#Model 0.11s var:269/1306 constraints:308/389
#Model 0.11s var:249/1306 constraints:291/389
#Model 0.12s var:244/1306 constraints:284/389
#Done 0.12s default_lp
#Done 0.12s quick_restart_no_lp
#Model 0.12s var:227/1306 constraints:261/389
#Model 0.12s var:213/1306 constraints:246/389
#Done 0.12s max_lp
Task timing n [ min, max] avg dev time n [ min, max] avg dev dtime
'core': 1 [ 70.68ms, 70.68ms] 70.68ms 0.00ns 70.68ms 1 [101.86ms, 101.86ms] 101.86ms 0.00ns 101.86ms
'default_lp': 1 [ 70.14ms, 70.14ms] 70.14ms 0.00ns 70.14ms 1 [ 43.82ms, 43.82ms] 43.82ms 0.00ns 43.82ms
'feasibility_pump': 3 [ 1.42ms, 10.03ms] 4.96ms 3.67ms 14.89ms 2 [ 1.51ms, 4.99ms] 3.25ms 1.74ms 6.50ms
'fj_long_default': 1 [ 17.44ms, 17.44ms] 17.44ms 0.00ns 17.44ms 1 [ 6.99ms, 6.99ms] 6.99ms 0.00ns 6.99ms
'fj_short_default': 1 [ 70.30ms, 70.30ms] 70.30ms 0.00ns 70.30ms 1 [ 1.54ms, 1.54ms] 1.54ms 0.00ns 1.54ms
'fs_random': 1 [ 16.62ms, 16.62ms] 16.62ms 0.00ns 16.62ms 1 [ 23.23us, 23.23us] 23.23us 0.00ns 23.23us
'graph_arc_lns': 2 [ 1.76ms, 11.16ms] 6.46ms 4.70ms 12.92ms 1 [351.78us, 351.78us] 351.78us 0.00ns 351.78us
'graph_cst_lns': 2 [ 6.52ms, 9.43ms] 7.97ms 1.45ms 15.95ms 2 [517.00ns, 1.49ms] 743.03us 742.52us 1.49ms
'graph_dec_lns': 2 [ 5.75ms, 15.99ms] 10.87ms 5.12ms 21.74ms 2 [489.00ns, 65.19us] 32.84us 32.35us 65.68us
'graph_var_lns': 2 [ 1.42ms, 2.74ms] 2.08ms 661.00us 4.16ms 1 [ 10.00ns, 10.00ns] 10.00ns 0.00ns 10.00ns
'max_lp': 1 [ 71.99ms, 71.99ms] 71.99ms 0.00ns 71.99ms 1 [ 5.54ms, 5.54ms] 5.54ms 0.00ns 5.54ms
'no_lp': 1 [ 70.28ms, 70.28ms] 70.28ms 0.00ns 70.28ms 1 [ 80.49ms, 80.49ms] 80.49ms 0.00ns 80.49ms
'pseudo_costs': 1 [ 70.89ms, 70.89ms] 70.89ms 0.00ns 70.89ms 1 [ 7.27ms, 7.27ms] 7.27ms 0.00ns 7.27ms
'quick_restart': 1 [ 70.30ms, 70.30ms] 70.30ms 0.00ns 70.30ms 1 [ 36.38ms, 36.38ms] 36.38ms 0.00ns 36.38ms
'quick_restart_no_lp': 1 [ 70.10ms, 70.10ms] 70.10ms 0.00ns 70.10ms 1 [ 75.31ms, 75.31ms] 75.31ms 0.00ns 75.31ms
'reduced_costs': 1 [ 71.34ms, 71.34ms] 71.34ms 0.00ns 71.34ms 1 [ 6.36ms, 6.36ms] 6.36ms 0.00ns 6.36ms
'rins/rens': 2 [ 1.32ms, 5.03ms] 3.17ms 1.86ms 6.35ms 1 [196.00us, 196.00us] 196.00us 0.00ns 196.00us
'rnd_cst_lns': 2 [ 3.79ms, 11.88ms] 7.83ms 4.04ms 15.66ms 2 [ 53.27us, 247.75us] 150.51us 97.24us 301.02us
'rnd_var_lns': 2 [ 1.71ms, 3.34ms] 2.53ms 813.00us 5.05ms 2 [ 3.09us, 41.30us] 22.19us 19.10us 44.39us
'routing_path_lns': 1 [ 7.44ms, 7.44ms] 7.44ms 0.00ns 7.44ms 1 [421.15us, 421.15us] 421.15us 0.00ns 421.15us
'routing_random_lns': 1 [ 11.39ms, 11.39ms] 11.39ms 0.00ns 11.39ms 1 [187.80us, 187.80us] 187.80us 0.00ns 187.80us
'violation_ls': 1 [ 49.23ms, 49.23ms] 49.23ms 0.00ns 49.23ms 1 [ 4.26ms, 4.26ms] 4.26ms 0.00ns 4.26ms
Search stats Bools Conflicts Branches Restarts BoolPropag IntegerPropag
'core': 1'269 2'631 18'331 4'354 567'096 502'679
'default_lp': 1'278 311 13'818 3'885 165'677 290'418
'fs_random': 1'255 84 2'422 2'422 92'129 81'230
'max_lp': 1'255 84 2'422 2'422 92'129 81'233
'no_lp': 1'255 2'071 32'869 5'061 392'175 497'149
'pseudo_costs': 1'255 84 2'422 2'422 93'109 83'178
'quick_restart': 1'300 172 13'949 3'851 150'297 262'154
'quick_restart_no_lp': 1'296 1'195 39'989 5'195 394'293 516'619
'reduced_costs': 1'255 84 2'422 2'422 93'107 83'200
SAT stats ClassicMinim LitRemoved LitLearned LitForgotten Subsumed MClauses MDecisions MLitTrue MSubsumed MLitRemoved MReused
'core': 2'153 31'058 42'280 0 52 718 9'933 34 462 4'274 42
'default_lp': 184 1'400 2'746 0 8 203 9'520 0 3 18 0
'fs_random': 0 0 84 0 0 0 0 0 0 0 0
'max_lp': 0 0 84 0 0 0 0 0 0 0 0
'no_lp': 1'673 17'082 37'496 0 24 212 23'779 0 8 99 4
'pseudo_costs': 0 0 84 0 0 0 0 0 0 0 0
'quick_restart': 62 326 894 0 0 196 8'865 0 0 0 1
'quick_restart_no_lp': 773 8'570 25'182 0 15 217 24'722 0 10 168 10
'reduced_costs': 0 0 84 0 0 0 0 0 0 0 0
Lp stats Component Iterations AddedCuts OPTIMAL DUAL_F. DUAL_U.
'default_lp': 3 824 24 345 0 2
'fs_random': 3 94 0 6 0 0
'max_lp': 1 271 401 9 0 0
'pseudo_costs': 1 371 446 19 0 0
'quick_restart': 3 520 23 262 0 4
'reduced_costs': 1 352 523 19 0 0
Lp dimension Final dimension of first component
'default_lp': 91 rows, 1238 columns, 242 entries with magnitude in [1.000000e+00, 1.000000e+00]
'fs_random': 103 rows, 1238 columns, 2337 entries with magnitude in [1.000000e+00, 1.000000e+00]
'max_lp': 627 rows, 1310 columns, 23337 entries with magnitude in [4.168404e-04, 1.000000e+00]
'pseudo_costs': 429 rows, 1310 columns, 5280 entries with magnitude in [2.436688e-05, 1.000000e+00]
'quick_restart': 95 rows, 1238 columns, 395 entries with magnitude in [1.000000e+00, 1.000000e+00]
'reduced_costs': 465 rows, 1310 columns, 5768 entries with magnitude in [1.902083e-04, 1.000000e+00]
Lp debug CutPropag CutEqPropag Adjust Overflow Bad BadScaling
'default_lp': 0 0 198 0 193 0
'fs_random': 0 0 2 0 0 0
'max_lp': 0 30 9 0 1'913 0
'pseudo_costs': 0 5 16 0 2'854 0
'quick_restart': 0 0 150 0 189 0
'reduced_costs': 0 1 16 0 2'500 0
Lp pool Constraints Updates Simplif Merged Shortened Split Strenghtened Cuts/Call
'default_lp': 177 90 755 100 639 0 0 24/98
'fs_random': 153 0 0 100 95 0 0 0/0
'max_lp': 823 48 4 100 100 0 385 401/716
'pseudo_costs': 868 200 991 100 964 5 140 446/913
'quick_restart': 176 124 926 100 745 0 0 23/100
'reduced_costs': 945 185 1'451 100 1'298 3 131 523/938
Lp Cut default_lp max_lp quick_restart reduced_costs pseudo_costs
CG_FF: 16 57 15 50 40
CG_K: - 26 - 24 -
CG_KL: - 1 - - -
CG_R: - 75 - 70 62
CG_RB: - 75 - 154 122
CG_RBP: - 9 - 18 5
Circuit: - 3 - 14 14
CircuitExact: - - - - 1
Clique: - - - 1 2
IB: - 27 - 57 59
MIR_1_FF: - 12 - 16 12
MIR_1_K: - 5 - 4 4
MIR_1_R: - 2 - 3 3
MIR_1_RB: - 9 - 18 7
MIR_1_RBP: - - - 2 2
MIR_2_FF: - 5 - 2 5
MIR_2_K: - 1 - 1 3
MIR_2_R: - 9 - 4 5
MIR_2_RB: - 5 - 3 3
MIR_2_RBP: - 3 - - 1
MIR_3_FF: - 2 - 4 4
MIR_3_K: - 1 - 1 1
MIR_3_R: - 7 - - 3
MIR_3_RB: - 2 - 4 4
MIR_3_RBP: - 6 - - 1
MIR_4_FF: - 1 - 2 5
MIR_4_K: - 1 - - -
MIR_4_R: - 1 - - 2
MIR_4_RB: - 3 - 2 1
MIR_4_RBP: - 2 - - -
MIR_5_FF: 3 7 6 11 10
MIR_5_K: - 1 - 1 1
MIR_5_R: - 3 - 1 1
MIR_5_RB: - 5 - 2 -
MIR_5_RBP: - 1 - - -
MIR_6_FF: - 3 - 5 4
MIR_6_R: - 2 - - 1
MIR_6_RB: - 1 - 3 2
ZERO_HALF_FF: 5 9 2 10 13
ZERO_HALF_K: - 1 - - -
ZERO_HALF_R: - 14 - 23 29
ZERO_HALF_RB: - 4 - 13 14
LNS stats Improv/Calls Closed Difficulty TimeLimit
'graph_arc_lns': 1/2 100% 0.81 0.10
'graph_cst_lns': 1/2 50% 0.54 0.10
'graph_dec_lns': 1/2 50% 0.54 0.10
'graph_var_lns': 1/2 100% 0.81 0.10
'rins/rens': 2/2 100% 0.81 0.10
'rnd_cst_lns': 0/2 100% 0.81 0.10
'rnd_var_lns': 0/2 100% 0.81 0.10
'routing_path_lns': 0/1 100% 0.71 0.10
'routing_random_lns': 0/1 100% 0.71 0.10
LS stats Batches Restarts LinMoves GenMoves CompoundMoves WeightUpdates
'fj_long_default': 1 1 9'511 0 0 489
'fj_short_default': 1 1 0 653 157 0
'violation_ls': 1 0 0 265 1 1
Solutions (2) Num Rank
'default_lp': 2 [1,2]
Objective bounds Num
'core': 1
'initial_domain': 1
'max_lp': 3
'pseudo_costs': 1
'quick_restart': 1
Solution repositories Added Queried Ignored Synchro
'feasible solutions': 14 29 0 11
'lp solutions': 0 0 0 0
'pump': 2 2
Improving bounds shared Num
'core': 143
'default_lp': 612
'no_lp': 18
'pseudo_costs': 2
'quick_restart': 271
'quick_restart_no_lp': 65
'reduced_costs': 9
Clauses shared Num
'core': 200
'default_lp': 10
'no_lp': 25
'quick_restart': 1
'quick_restart_no_lp': 102
CpSolverResponse summary:
status: OPTIMAL
objective: 551
best_bound: 551
integers: 1348
booleans: 1255
conflicts: 84
branches: 2422
propagations: 92129
integer_propagations: 81230
restarts: 2422
lp_iterations: 94
walltime: 0.135173
usertime: 0.135173
deterministic_time: 0.421991
gap_integral: 0.322107
solution_fingerprint: 0x9d0e69acfb4c19b3
status OPTIMALsol = [i for (v, i) in sorted([(solver.Value(u[i]), i) for i in u])]
if status == cp_model.OPTIMAL:
print('Optimal Value:', solver.ObjectiveValue())Optimal Value: 551.0import networkx as nx
G = nx.Graph()
for idx, i in enumerate(sol[:-1]):
G.add_edge(i, sol[idx + 1])
G.add_edge(sol[-1],sol[1])
pos = {i: (x_[i-1], y_[i-1]) for i in range(1,n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");