グラフ分割問題に対する定式化とメタヒューリスティクス

準備

ここでは,付録2で準備したグラフに関する基本操作を集めたモジュール graphtools.py を読み込んでいる. 環境によって,モジュールファイルの場所は適宜変更されたい.

グラフ2分割問題

グラフ2分割問題(graph bipartition problem)は、以下のように定義される。

点数 $n=|V|$ が偶数である無向グラフ $G=(V,E)$ が与えられたとき, 点集合 $V$ の等分割(uniform partition, eqipartition) $(L,R)$ とは, $L \cap R = \emptyset$, $L \cup R =V$, $|L|=|R|=n/2$ を満たす点の部分集合の対である.

$L$ と $R$ をまたいでいる枝 (正確には $i \in L, j \in R$ もしくは $i \in R, j \in L$ を満たす枝 $(i,j)$) の本数を最小にする等分割 $(L,R)$ を求める問題がグラフ2分割問題である.

一般には $k (>2)$ 分割も考えられるが,ここでは $k=2$ に絞って議論する.

問題を明確化するために,グラフ分割問題を整数最適化問題として定式化しておく. 無向グラフ $G=(V,E)$ に対して,$L \cap R =\emptyset$(共通部分がない), $L \cup R =V$(合わせると点集合全体になる)を満たす非順序対 $(L,R)$ を2分割(bipartition)と呼ぶ. 分割 $(L,R)$ において,$L$ は左側,$R$ は右側を表すが,これらは逆にしても同じ分割であるので,非順序対と呼ばれる.

点 $i$ が,分割 $(L,R)$ の $L$ 側に含まれているとき $1$, それ以外の($R$ 側に含まれている)とき $0$ の $0$-$1$ 変数 $x_i$ を導入する. このとき,等分割であるためには,$x_i$ の合計が $n/2$ である必要がある. 枝 $(i,j)$ が $L$ と $R$ をまたいでいるときには,$x_i (1-x_j)$ もしくは $(1-x_i)x_j$ が $1$ になることから,以下の定式化を得る.

$$ \begin{array}{l l l} minimize & \sum_{(i,j) \in E} x_i (1-x_j)+(1-x_i)x_j & \\ s.t. & \sum_{i \in V} x_i = n/2 & \\ & x_i \in \{0,1\} & \forall i \in V \end{array} $$

数理最適化ソルバーの多くは,上のように凸でない2次の項を目的関数に含んだ最小化問題には対応していない. したがって,一般的な(混合)整数最適化ソルバーで求解するためには,2次の項を線形関数に変形してから解く必要がある.

枝 $(i,j)$ が $L$ と $R$ をまたいでいるとき $1$,それ以外のとき $0$ になる $0$-$1$ 変数 $y_{ij}$ を導入する. すると,上の2次整数最適化問題は,以下の線形整数最適化に帰着される.

$$ \begin{array}{l l l} minimize & \sum_{(i,j) \in E} y_{ij} & \\ s.t. & \sum_{i \in V} x_i = n/2 & \\ & x_i -x_j \leq y_{ij} & \forall (i,j) \in E \\ & x_j -x_i \leq y_{ij} & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \\ & y_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array} $$

最初の制約は等分割を規定する. 2番目の制約は,$i \in L$ で $j \not\in L$ のとき $y_{ij}=1$ になることを規定する. 3番目の制約は,$j \in L$ で $i \not\in L$ のとき $y_{ij}=1$ になることを規定する.

上の定式化を行う関数を以下に示す.

Pythonによる求解

make_data[source]

make_data(n, prob)

make_data: prepare data for a random graph Parameters:

- n: number of vertices
- prob: probability of existence of an edge, for each pair of vertices

Returns a tuple with a list of vertices and a list edges.

gpp[source]

gpp(V, E)

gpp -- model for the graph partitioning problem Parameters:

- V: set/list of nodes in the graph
- E: set/list of edges in the graph

Returns a model, ready to be solved.

V, E = make_data(10, 0.4)
model = gpp(V, E)
model.optimize()
print("Opt.value=", model.ObjVal)
x, y = model.__data
L = [i for i in V if x[i].X >= 0.5]
R = [i for i in V if x[i].X < 0.5]
Set parameter Username
Academic license - for non-commercial use only - expires 2023-11-07
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 45 rows, 32 columns and 142 nonzeros
Model fingerprint: 0x64d2af4c
Variable types: 0 continuous, 32 integer (32 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [5e+00, 5e+00]
Found heuristic solution: objective 15.0000000
Presolve time: 0.00s
Presolved: 45 rows, 32 columns, 142 nonzeros
Variable types: 0 continuous, 32 integer (32 binary)

Root relaxation: objective 0.000000e+00, 15 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    0.00000    0   10   15.00000    0.00000   100%     -    0s
H    0     0                      11.0000000    0.00000   100%     -    0s
     0     0    5.00000    0   10   11.00000    5.00000  54.5%     -    0s
H    0     0                      10.0000000    5.00000  50.0%     -    0s
H    0     0                       9.0000000    5.00000  44.4%     -    0s
H    0     0                       8.0000000    5.00000  37.5%     -    0s

Cutting planes:
  Gomory: 3
  MIR: 2
  Zero half: 4
  RLT: 5

Explored 1 nodes (43 simplex iterations) in 0.02 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)

Solution count 5: 8 9 10 ... 15

Optimal solution found (tolerance 1.00e-04)
Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0000%
Opt.value= 8.0
G = nx.Graph()
G.add_edges_from(E)
pos = nx.layout.spring_layout(G)
nx.draw(G, pos=pos, with_labels=True, node_size=500, node_color="yellow", nodelist=L)
edgelist = [(i, j) for (i, j) in E if (i in L and j in R) or (i in R and j in L)]
nx.draw(
    G,
    pos=pos,
    with_labels=False,
    node_size=500,
    nodelist=R,
    node_color="Orange",
    edgelist=edgelist,
    edge_color="red",
)

AMPLによる求解

AMPLモデル(gpp.mod)

set V;
set E within {V,V};
param n := card(V);

var x {V} binary;
var y {E} binary;

minimize z: sum {(i,j) in E} y[i,j];

subject to
Equipart: sum {i in V} x[i] = n/2;
AcrossLR {(i,j) in E}: x[i] - x[j] <= y[i,j];
AcrossRL {(i,j) in E}: x[j] - x[i] <= y[i,j];
ampl = AMPL(Environment(AMPL_PATH))
ampl.option['solver'] = 'highs'
ampl.read("../ampl/gpp.mod")
ampl.set['V'] = list(V)
ampl.set['E'] = list(E)

ampl.solve()
z = ampl.obj['z']
print("Optimum:", z.value())

x = ampl.var['x']
L = [i for i in V if x[i].value() >.5]
R = [i for i in V if x[i].value() <.5]
print("L =", L)
print("R =", R)

y = ampl.var['y']
across = [(i,j) for (i,j) in E if y[i,j].value() > .5]
print("Across edges:", across)
HiGHS 1.2.2:HiGHS 1.2.2: optimal solution; objective 8
3 branching nodes
Optimum: 8.0
L = [1, 3, 5, 9, 10]
R = [2, 4, 6, 7, 8]
Across edges: [(1, 2), (1, 6), (1, 7), (1, 8), (2, 3), (2, 5), (5, 7), (8, 10)]

タブーサーチ

タブーサーチ(tabu search)は,Gloverによって提案されたメタヒューリスティクスであり, グラフ分割問題に対しては,点をLからRに,RからLに移す簡単な近傍をもとに設計できる. 一度移動させた点は,しばらく移動させないことによって,同じ解に戻らないようにすることが,タブー(禁断)の名前の由来である.

以下のメタヒューリスティクスの詳細については,拙著「メタヒューリスティクスの数理」(共立出版)を参照されたい.

construct[source]

construct(nodes)

A simple construction method.

The solution is represented by a vector:

  • sol[i]=0 - node i is in partition 0
  • sol[i]=1 - node i is in partition 1

evaluate[source]

evaluate(nodes, adj, sol, alpha=0.0)

Evaluate a solution.

Input:

  • nodes, adj - the instance's graph
  • sol - the solution to evaluate
  • alpha - a penalty for imbalanced solutions

Determines:

  • the cost of a solution, i.e., the number of edges going from one partition to the other;
  • bal - balance, i.e., the number of vertices in excess in partition 0
  • s[i] - number of edges adjacent to i in the same partition;
  • d[i] - number of edges adjacent to i in a different partition.

find_move_rnd[source]

find_move_rnd(part, nodes, adj, sol, s, d, tabu, tabulen, iteration)

Probabilistically find the best non-tabu move into partition type 'part'.

find_move[source]

find_move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration)

Find the best non-tabu move into partition type 'part'.

move[source]

move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration)

Determine and execute the best non-tabu move.

tabu_search(nodes, adj, sol, max_iter, tabulen, report=None)

Execute a tabu search run.

rndseed = 2
random.seed(rndseed)
n = 100
pos = {i: (random.random(), random.random()) for i in range(n)}
G = nx.random_geometric_graph(n, 0.2, pos=pos)

nodes = G.nodes()
adj = [set([]) for i in nodes]
for (i, j) in G.edges():
    adj[i].add(j)
    adj[j].add(i)

max_iterations = 1000
tabulen = len(nodes) / 2

sol = construct(nodes)
z, _, s, d = evaluate(nodes, adj, sol)
print("initial partition: z =", z)
print(sol)
print()
print("starting tabu search,", max_iterations, "iterations, tabu length =", tabulen)
sol, cost = tabu_search(nodes, adj, sol, max_iterations, tabulen)
print("final solution: z=", cost)
print(sol)
initial partition: z = 239.0
[0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1]

starting tabu search, 1000 iterations, tabu length = 50.0
final solution: z= 31.0
[1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1]
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
    (i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
    G,
    pos=pos,
    with_labels=False,
    node_size=100,
    nodelist=R,
    node_color="Orange",
    edgelist=edgelist,
    edge_color="red",
)

アニーリング法

ランダム性を用いて局所的最適解からの脱出を試みるメタヒューリスティクスの代表例として アニーリング法(simulated annealing method)がある. この解法の土台は,1953年に Metropolis--Rosenbluth--Rosenbluth--Teller--Tellerによって築かれたものであり,何度も再発見されている. そのため,様ざなま別名をもつ.一部を紹介すると, モンテカルロ焼なまし法(Monte Carlo annealing), 確率的丘登り法(probabilistic hill climbing), 統計的冷却法(statistical cooling), 確率的緩和法(stochastic relaxation), Metropolisアルゴリズム(Metropolis algorithm) などがある.

アニーリング法の特徴は,温度とよばれるパラメータを徐々に小さくしていくことによって,改悪の確率を $0$ に近づけていくことである.. 以下では,分割の非均衡度 $|L|-|R|$ をペナルティとしたアニーリング法を示す.

find_move_rnd_sa[source]

find_move_rnd_sa(n, sol, alpha, s, d, bal)

Find a random node to move from one part into the other. For simulated annealing.

update_move[source]

update_move(adj, sol, s, d, bal, istar)

Execute the chosen move.

metropolis[source]

metropolis(T, delta)

Metropolis criterion for new configuration acceptance

estimate_temperature[source]

estimate_temperature(n, sol, s, d, bal, X0, alpha)

Estimate initial temperature: check empirically based on a series of 'ntrials', that the estimated temperature leads to a rate 'X0'% acceptance:

annealing[source]

annealing(nodes, adj, sol, initprob, L, tempfactor, freezelim, minpercent, alpha, report)

Simulated annealing for the graph partitioning problem

Parameters:

  • nodes, adj - graph definition
  • sol - initial solution
  • initprob - initial acceptance rate
  • L - number of tentatives at each temperature
  • tempfactor - cooling ratio
  • freezelim - max number of iterations with less that minpercent acceptances
  • minpercent - percentage of accepted moves for being not frozen
  • report - function used for output of best found solutions
initprob = 0.4  # initial acceptance probability
sizefactor = 16  # for det. # tentatives at each temp.
tempfactor = 0.95  # cooling ratio
freezelim = 5  # max number of iterations with less that minpercent acceptances
minpercent = 0.02  # fraction of accepted moves for being not frozen
alpha = 3.0  # imballance factor
N = len(nodes)  # neighborhood size
L = N * sizefactor  # number of tentatives at current temperature

print("starting simulated annealing, parameters:")
print("   initial acceptance probability", initprob)
print("   cooling ratio", tempfactor)
print("   # tentatives at each temp.", L)
print("   percentage of accepted moves for being not frozen", minpercent)
print("   max # of it.with less that minpercent acceptances", freezelim)
print("   imballance factor", alpha)
print()
sol = construct(nodes)  # starting solution
z, bal, s, d = evaluate(nodes, adj, sol, alpha)
print("initial partition: z =", z)
print(sol)
print()

LOG = False
sol, z = annealing(
    nodes, adj, sol, initprob, L, tempfactor, freezelim, minpercent, alpha, report=False
)
zp, bal, sp, dp = evaluate(nodes, adj, sol, alpha)
assert abs(zp - z) < 1.0e-9  # floating point approx.equality

print()
print("final solution: z=", z)
print(sol)
starting simulated annealing, parameters:
   initial acceptance probability 0.4
   cooling ratio 0.95
   # tentatives at each temp. 1600
   percentage of accepted moves for being not frozen 0.02
   max # of it.with less that minpercent acceptances 5
   imballance factor 3.0

initial partition: z = 238.0
[0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1]


final solution: z= 33.0
[0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0]
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
    (i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
    G,
    pos=pos,
    with_labels=False,
    node_size=100,
    nodelist=R,
    node_color="Orange",
    edgelist=edgelist,
    edge_color="red",
)

集中化と多様化を入れたタブーサーチ

良いメタヒューリスティクスを設計するためのコツは,集中化と多様化の両者をバランス良く組み込むことである. ここで集中化(intensification)とは,良い解のまわりを集中して探索することであり, 多様化(diversification)とは,まだ探していない解を探索することである.

集中化と多様化を加味したタブーサーチを以下に示す.

diversify[source]

diversify(sol, nodes)

Diversify: keep a part of the solution with random size.

The solution is represented by a vector:

  • sol[i]=0 - node i is in partition 0
  • sol[i]=1 - node i is in partition 1

ts_intens_divers[source]

ts_intens_divers(nodes, adj, sol, max_iter, tabulen, report)

Execute a tabu search run, with intensification/diversification.

sol = construct(nodes)
print("initial partition: z =", z)
print(sol)
print()
print("starting tabu search,", max_iterations, "iterations, tabu length =", tabulen)
sol, cost = ts_intens_divers(nodes, adj, sol, max_iterations, tabulen, report=False)
print("final solution: z=", cost)
print(sol)
initial partition: z = 33.0
[0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0]

starting tabu search, 1000 iterations, tabu length = 50.0
final solution: z= 31.0
[1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0]
#     G.nodes[v]["color"] = sol[v]
# fig = gts.to_plotly_fig(G,node_size=30,
#                     line_width=0.1,
#                     line_color='blue',
#                     text_size=10, pos=pos)
# plotly.offline.plot(fig);
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
    (i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
    G,
    pos=pos,
    with_labels=False,
    node_size=100,
    nodelist=R,
    node_color="Orange",
    edgelist=edgelist,
    edge_color="red",
)

グラフ多分割問題

上では,グラフの2分割問題を考えたが,3以上の複数の集合に,できるだけ均等に分割することを考える. この問題は,グラフ多分割問題(multiway graph partitioning problem)とよばれ,様々な制約が付加された問題が考えられている.

この問題に対するメタヒューリスティクスが,以下のプロジェクトで管理されている.

このパッケージを使うことによって,様々な付加制約がついたグラフ多分割問題を,高速に解くことができる.

最大カット問題

最大カット問題は、以下のように定義される。

無向グラフ $G=(V,E)$ が与えられたとき, 点集合 $V$ の部分集合 $S$ で, 両端点が $S $ と $V \setminus S$ に含まれる枝集合(カット)の位数が最大のものをもとめよ.

もしくは,枝 $(i,j) \in E$ の上に重み $w_{ij}$ が定義されていて,カットに含まれる枝の重みの合計を最大化する.

メタヒューリスティクスは,グラフ分割問題と同様に設計できる.

線形定式化

グラフ分割問題と同様に, 点 $i$ が,分割 $(L,R)$ の $L$ 側に含まれているとき $1$, それ以外の($R$ 側に含まれている)とき $0$ の $0$-$1$ 変数 $x_i$ を導入する. 枝 $(i,j)$ が $L$ と $R$ をまたいでいるとき $1$,それ以外のとき $0$ になる $0$-$1$ 変数 $y_{ij}$ を導入する.

2分割を表す制約のかわりに, 点 $i,j$ が同じ側に含まれているときに $y_{ij}$ を $0$ にするための制約を追加する.

$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} w_{ij} y_{ij} & \\ s.t. & x_i + x_j \geq y_{ij} & \forall (i,j) \in E \\ & 2- x_i - x_j \geq y_{ij} & \forall (i,j) \in E \\ & x_i -x_j \leq y_{ij} & \forall (i,j) \in E \\ & x_j -x_i \leq y_{ij} & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \\ & y_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array} $$
def maxcut(V, E):
    """maxcut -- model for the graph maxcut problem
    Parameters:
        - V: set/list of nodes in the graph
        - E: set/list of edges in the graph
    Returns a model, ready to be solved.
    """
    model = Model("maxcut")
    x = {}
    y = {}
    for i in V:
        x[i] = model.addVar(vtype="B", name=f"x({i})")
    for (i, j) in E:
        y[i, j] = model.addVar(vtype="B", name=f"y({i},{j})")
    model.update()

    for (i, j) in E:
        model.addConstr(x[i] + x[j] >= y[i, j], f"Edge({i},{j})")
        model.addConstr(2 - x[j] - x[i] >= y[i, j], f"Edge({j},{i})")
        model.addConstr(x[i] - x[j] <= y[i, j], f"EdgeLB1({i},{j})")
        model.addConstr(x[j] - x[i] <= y[i, j], f"EdgeLB2({j},{i})")

    model.setObjective(quicksum(y[i, j] for (i, j) in E), GRB.MAXIMIZE)

    model.update()
    model.__data = x
    return model

2次錐最適化による定式化

ここでは,以下の論文による2次錐最適化による強い定式化を示す.

A NEW SECOND-ORDER CONE PROGRAMMING RELAXATION FOR MAX-CUT PROBLEMS, Masakazu Muramatsu Tsunehiro Suzuki, Journal of the Operations Research, 2003, Vol. 46, No. 2, 164-177

$x_j$ は上の定式化と同じであるが,以下の2つの $0$-$1$ 変数を導入する.

  • $s_{ij}$: $i$ と $j$ が同じ分割に含まれるとき $1$
  • $z_{ij}$: $i$ と $j$ が異なる分割に含まれるとき $1$
$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} w_{ij} z_{ij} & \\ s.t. & (x_i + x_j-1)^2 \leq s_{ij} & \forall (i,j) \in E \\ & (x_j- x_i)^2 \leq z_{ij} & \forall (i,j) \in E \\ & s_{ij} + z_{ij} = 1 & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \\ & s_{ij} \in \{0,1\} & \forall (i,j) \in E \\ & z_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array} $$
def maxcut_soco(V, E):
    """maxcut_soco -- model for the graph maxcut problem
    Parameters:
        - V: set/list of nodes in the graph
        - E: set/list of edges in the graph
    Returns a model, ready to be solved.
    """
    model = Model("max cut -- scop")
    x, s, z = {}, {}, {}
    for i in V:
        x[i] = model.addVar(vtype="B", name=f"x({i})")
    for (i, j) in E:
        s[i, j] = model.addVar(vtype="C", name=f"s({i},{j})")
        z[i, j] = model.addVar(vtype="C", name=f"z({i},{j})")
    #model.update()

    for (i, j) in E:
        model.addConstr((x[i] + x[j] - 1) * (x[i] + x[j] - 1) <= s[i, j], f"S({i},{j})")
        model.addConstr((x[j] - x[i]) * (x[j] - x[i]) <= z[i, j], f"Z({i},{j})")
        model.addConstr(s[i, j] + z[i, j] == 1, f"P({i},{j})")

    model.setObjective(quicksum(z[i, j] for (i, j) in E), GRB.MAXIMIZE)

    #model.update()
    #model.__data = x
    return model

def make_data(n, prob):
    """make_data: prepare data for a random graph
    Parameters:
        - n: number of vertices
        - prob: probability of existence of an edge, for each pair of vertices
    Returns a tuple with a list of vertices and a list edges.
    """
    V = range(1, n + 1)
    E = [(i, j) for i in V for j in V if i < j and random.random() < prob]
    return V, E
#GRB = MDO
from gurobipy import Model, quicksum, GRB

random.seed(3)
V, E = make_data(100, 0.1)

model = maxcut(V, E)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
    print("Opt.value=", model.ObjVal)
    obj_lin = model.ObjVal
    x = model.__data
    print([i for i in V if x[i].X >= 0.5])
    print([i for i in V if x[i].X < 0.5])
else:
    pass
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 1948 rows, 587 columns and 5844 nonzeros
Model fingerprint: 0xa1aec4f7
Variable types: 0 continuous, 587 integer (587 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+00, 2e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 974 rows and 0 columns
Presolve time: 0.01s
Presolved: 974 rows, 587 columns, 2922 nonzeros
Variable types: 0 continuous, 587 integer (587 binary)
Found heuristic solution: objective 274.0000000

Root relaxation: objective 4.870000e+02, 108 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  487.00000    0  100  274.00000  487.00000  77.7%     -    0s
H    0     0                     309.0000000  487.00000  57.6%     -    0s
H    0     0                     314.0000000  394.64321  25.7%     -    0s
H    0     0                     317.0000000  394.64321  24.5%     -    0s
H    0     0                     318.0000000  394.64321  24.1%     -    0s
     0     0  394.64321    0  337  318.00000  394.64321  24.1%     -    0s
H    0     0                     320.0000000  380.34427  18.9%     -    0s
H    0     0                     325.0000000  380.34427  17.0%     -    0s
     0     0  380.34427    0  452  325.00000  380.34427  17.0%     -    0s
H    0     0                     326.0000000  378.28358  16.0%     -    0s
H    0     0                     329.0000000  378.28358  15.0%     -    0s
     0     0  378.28358    0  454  329.00000  378.28358  15.0%     -    0s
     0     0  377.59775    0  470  329.00000  377.59775  14.8%     -    0s
     0     0  377.46366    0  477  329.00000  377.46366  14.7%     -    0s
     0     0  377.41514    0  483  329.00000  377.41514  14.7%     -    0s
     0     0  374.40557    0  493  329.00000  374.40557  13.8%     -    0s
     0     0  373.83797    0  509  329.00000  373.83797  13.6%     -    0s
     0     0  373.68816    0  515  329.00000  373.68816  13.6%     -    0s
     0     0  373.60416    0  523  329.00000  373.60416  13.6%     -    0s
H    0     0                     330.0000000  372.51558  12.9%     -    1s
H    0     0                     331.0000000  372.51558  12.5%     -    1s
     0     0  372.51558    0  526  331.00000  372.51558  12.5%     -    1s
H    0     0                     346.0000000  372.51558  7.66%     -    1s
     0     0  372.16728    0  533  346.00000  372.16728  7.56%     -    1s
     0     0  372.04505    0  532  346.00000  372.04505  7.53%     -    1s
     0     0  371.99409    0  535  346.00000  371.99409  7.51%     -    1s
     0     0  370.83131    0  525  346.00000  370.83131  7.18%     -    1s
     0     0  370.32761    0  535  346.00000  370.32761  7.03%     -    1s
     0     0  370.08554    0  535  346.00000  370.08554  6.96%     -    1s
     0     0  369.94819    0  551  346.00000  369.94819  6.92%     -    1s
     0     0  369.89772    0  556  346.00000  369.89772  6.91%     -    1s
     0     0  368.95515    0  544  346.00000  368.95515  6.63%     -    2s
H    0     0                     348.0000000  368.63515  5.93%     -    2s
H    0     0                     349.0000000  368.63515  5.63%     -    2s
     0     0  368.63515    0  552  349.00000  368.63515  5.63%     -    2s
     0     0  368.44040    0  554  349.00000  368.44040  5.57%     -    2s
     0     0  368.37365    0  556  349.00000  368.37365  5.55%     -    2s
     0     0  367.59563    0  549  349.00000  367.59563  5.33%     -    2s
     0     0  367.30850    0  557  349.00000  367.30850  5.25%     -    2s
     0     0  367.17035    0  562  349.00000  367.17035  5.21%     -    2s
     0     0  367.06203    0  573  349.00000  367.06203  5.18%     -    2s
     0     0  366.72354    0  564  349.00000  366.72354  5.08%     -    2s
     0     0  366.72354    0  562  349.00000  366.72354  5.08%     -    2s
     0     2  366.72354    0  562  349.00000  366.72354  5.08%     -    3s
    63    72  364.04029    7  562  349.00000  365.52425  4.73%  1173    5s
H  515   503                     350.0000000  365.52425  4.44%   576    8s
   652   622  364.62803    7  562  350.00000  365.27255  4.36%   541   10s
H  689   531                     353.0000000  365.27255  3.48%   538   10s
H  700   506                     354.0000000  365.27255  3.18%   539   10s
H 1017   711                     355.0000000  365.14244  2.86%   528   13s
  1178   847  358.44811   26  520  355.00000  365.14244  2.86%   540   15s
  1494  1002  363.54941    9  680  355.00000  364.98151  2.81%   536   20s
  1675  1141  359.50612   17  677  355.00000  362.43523  2.09%   567   25s
  2003  1214  356.67364   25  659  355.00000  362.43523  2.09%   557   30s
  2533  1280  357.77001   20  666  355.00000  361.90278  1.94%   536   35s
  3163  1344  357.00634   23  668  355.00000  361.18835  1.74%   515   40s

Cutting planes:
  MIR: 17
  Flow cover: 13

Explored 3785 nodes (1938869 simplex iterations) in 44.84 seconds (92.45 work units)
Thread count was 16 (of 16 available processors)

Solution count 10: 355 354 353 ... 329

Solve interrupted
Best objective 3.550000000000e+02, best bound 3.600000000000e+02, gap 1.4085%
#GRB = MDO
from gurobipy import Model, quicksum, GRB

model = maxcut_soco(V, E)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
    print("Opt.value=", model.ObjVal)
    x = model.__data
    print([i for i in V if x[i].X >= 0.5])
    print([i for i in V if x[i].X < 0.5])
    # assert model.ObjVal == obj_lin
else:
    pass
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 487 rows, 1074 columns and 974 nonzeros
Model fingerprint: 0x9f4edd72
Model has 974 quadratic constraints
Variable types: 974 continuous, 100 integer (100 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 2e+00]
  QLMatrix range   [1e+00, 2e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
  QRHS range       [1e+00, 1e+00]
Presolve time: 0.01s
Presolved: 2922 rows, 1561 columns, 8279 nonzeros
Variable types: 974 continuous, 587 integer (587 binary)

Interrupt request received
Found heuristic solution: objective -0.0000000

Root relaxation: objective 4.870000e+02, 118 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  487.00000    0  100   -0.00000  487.00000      -     -    0s
H    0     0                     286.0000000  487.00000  70.3%     -    0s
H    0     0                     289.0000000  487.00000  68.5%     -    0s
H    0     0                     297.0000000  487.00000  64.0%     -    0s
H    0     0                     305.0000000  478.00000  56.7%     -    0s
     0     0  414.52500    0  290  305.00000  414.52500  35.9%     -    0s
     0     0  413.84615    0  283  305.00000  413.84615  35.7%     -    0s
H    0     0                     307.0000000  413.84615  34.8%     -    0s
     0     0  408.00000    0  179  307.00000  408.00000  32.9%     -    0s
     0     0  407.50000    0  184  307.00000  407.50000  32.7%     -    0s
H    0     0                     308.0000000  407.50000  32.3%     -    0s
     0     0  403.92188    0  366  308.00000  403.92188  31.1%     -    0s
     0     0  403.56897    0  339  308.00000  403.56897  31.0%     -    0s
H    0     0                     310.0000000  403.56897  30.2%     -    0s
     0     0  402.42134    0  433  310.00000  402.42134  29.8%     -    0s
     0     0  402.42134    0  383  310.00000  402.42134  29.8%     -    0s
H    0     0                     311.0000000  402.25597  29.3%     -    0s
H    0     2                     312.0000000  402.25597  28.9%     -    0s
     0     2  402.25597    0  378  312.00000  402.25597  28.9%     -    0s
H   33    40                     317.0000000  393.87912  24.3%   374    0s
H   66    76                     320.0000000  393.87912  23.1%   312    1s
H   68    76                     322.0000000  393.87912  22.3%   312    1s
H  105   116                     323.0000000  393.87912  21.9%   271    1s
H  106   116                     324.0000000  393.87912  21.6%   270    1s
H  109   116                     325.0000000  393.87912  21.2%   267    1s
H  112   116                     326.0000000  393.87912  20.8%   265    1s
H  142   154                     330.0000000  393.87912  19.4%   245    1s
*  161   164              30     334.0000000  393.87912  17.9%   239    1s
*  179   165              28     336.0000000  393.87912  17.2%   226    1s
H  184   174                     337.0000000  393.87912  16.9%   222    1s
H  184   174                     340.0000000  393.87912  15.8%   222    1s
H  192   174                     343.0000000  393.87912  14.8%   220    1s
H  222   213                     347.0000000  393.87912  13.5%   221    1s
H  255   247                     353.0000000  393.87912  11.6%   210    1s
H  260   247                     354.0000000  393.87912  11.3%   208    1s
H  309   251                     355.0000000  391.69994  10.3%   190    1s

2次無制約2値最適化による定式化

Gurobiの最近のバージョンでは,2次無制約2値最適化 (Quadratic unconstrained binary optimization: QUBO)に対する専用のヒューリスティクスを組み込んでいる.

最大カット問題は,QUBOとして以下のように定式化できる. QUBOでは,2値変数は $-1$ か $1$ の値ととるものと仮定する. 分割 $(L,R)$ の $L$ 側に含まれているとき $1$, それ以外の($R$ 側に含まれている)とき $-1$ の2値変数 $x_i$ を用いる.

$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} w_{ij} (1-x_{i}x_{j}) & \\ & x_i \in \{-1,1\} & \forall i \in V \end{array} $$

これを,以下の変換を用いて $0$-$1$ 変数 $\hat{x}_i$ に変換する.

$$ \hat{x}_i = \frac{x_i+1}{2} $$

これをQUBOによる代入すると,2次無制約の$0$-$1$ 整数最適化は,以下のようになる.

$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} \frac{1}{2} w_{ij} (-4 \hat{x}_{i} \hat{x}_{j}+ 2\hat{x}_{i} + 2 \hat{x}_{j}) & \\ & \hat{x}_i \in \{0,1\} & \forall i \in V \end{array} $$

以下のサイトからベンチマーク問題例を入手できるので, ベンチマーク問題例で実験する.

http://grafo.etsii.urjc.es/optsicom/maxcut/

from gurobipy import Model, quicksum, GRB

folder = "../data/maxcut/"
gurobi = []
for iter in range(1,2): #ここで問題番号を指定する(全部で54問)
    f = open(folder + f"g{iter}.rud")
    data = f.readlines()
    f.close()

    n, m = list(map(int, data[0].split()))
    G = nx.Graph()
    edges = []
    for row in data[1:]:
        i, j, w = list(map(int, row.split()))
        edges.append((i, j, w))
    G.add_weighted_edges_from(edges)

    model = Model()

    x = {}
    for i in G.nodes():
        x[i] = model.addVar(name=f"x[{i}]",vtype="B")
    model.setObjective( quicksum(1/2* G[i][j]["weight"]*(-4*x[i]*x[j]+2*x[i]+2*x[j])  for (i, j) in G.edges()), GRB.MAXIMIZE)
    model.Params.TimeLimit = 10
    model.Params.OutputFlag=1
    model.optimize()
    print(iter, model.ObjVal)
    gurobi.append(model.ObjVal)
Set parameter TimeLimit to value 10
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 0 rows, 800 columns and 0 nonzeros
Model fingerprint: 0xe5eaeb92
Model has 19176 quadratic objective terms
Variable types: 0 continuous, 800 integer (800 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [3e+01, 7e+01]
  QObjective range [4e+00, 4e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]
Found heuristic solution: objective -0.0000000
Found heuristic solution: objective 11557.000000
Presolve time: 0.44s
Presolved: 0 rows, 800 columns, 0 nonzeros
Presolved model has 19976 quadratic objective terms
Variable types: 0 continuous, 800 integer (800 binary)
Found heuristic solution: objective 11578.000000

Root relaxation: objective 1.223167e+04, 1192 iterations, 0.75 seconds (0.73 work units)

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

     0     0 12231.6655    0  799 11578.0000 12231.6655  5.65%     -    1s
     0     0 12231.6655    0  799 11578.0000 12231.6655  5.65%     -    1s
     0     2 12231.6655    0  799 11578.0000 12231.6655  5.65%     -    2s
H    4     8                    11580.000000 12231.3774  5.63%   5.2    2s
H   29    32                    11599.000000 12230.9681  5.45%   3.4    2s
H   42    48                    11624.000000 12230.9681  5.22%   3.4    2s
   148   157 12221.7435   18  776 11624.0000 12228.9709  5.20%   3.8    5s
   399   424 12215.7193   40  739 11624.0000 12226.6652  5.18%   4.0   10s

Explored 423 nodes (2877 simplex iterations) in 10.02 seconds (6.70 work units)
Thread count was 16 (of 16 available processors)

Solution count 6: 11624 11599 11580 ... -0

Time limit reached
Best objective 1.162400000000e+04, best bound 1.222666515917e+04, gap 5.1847%
1 11624.0

制約最適化ソルバー SCOP による求解

上と同じ問題例に対して,制約最適化ソルバー SCOP (付録1参照)で求解を試みる.

変数は $X_i (i \in V)$ で値集合は $\{0,1\}$ とする. 値変数は $x_{i0}, x_{i1}$ となる.

ベンチマーク問題例には,枝の重み $w_{ij}, (i,j) \in E$が $1$ のものと $-1$ のものがある.

枝の重みが $1$ の枝に対しては, 重み(逸脱ペナルティ)が $1$ の考慮制約として,以下の2次制約を定義する.

$$ x_{i0} x_{j1} + x_{i1} x_{j0} \geq 1 \ \ \ \forall (i,j) \in E, w_{ij}=1 $$

枝 $(i,j)$ の両端点が異なる集合に含まれるとき左辺は $1$ となり,制約は満たされる. それ以外のとき左辺は $0$ となり,制約は $1$ だけ破られる.

枝の重みが $-1$ の枝に対しては, 重み(逸脱ペナルティ)が $1$ の考慮制約として,以下の2次制約を定義する.

$$ x_{i0} x_{j1} + x_{i1} x_{j0} \leq 0 \ \ \ \forall (i,j) \in E, w_{ij}=-1 $$

枝 $(i,j)$ の両端点が異なる集合に含まれるとき左辺は $1$ となり,制約は $1$ だけ破られる.

枝の本数 $|E|$ から逸脱の総数を減じたものが,目的関数になる.

SCOPモジュールからインポートする.

from scop import *
folder = "../data/maxcut/"
for iter in range(1, 2): #ここで問題番号を指定する(全部で54問)
    f = open(folder + f"g{iter}.rud")
    data = f.readlines()
    f.close()

    n, m = list(map(int, data[0].split()))
    G = nx.Graph()
    edges = []
    for row in data[1:]:
        i, j, w = list(map(int, row.split()))
        edges.append((i, j, w))
    G.add_weighted_edges_from(edges)

    model = Model()

    x = {}
    for i in G.nodes():
        x[i] = model.addVariable(name=f"x[{i}]", domain=[0, 1])
    for (i, j) in G.edges():
        if G[i][j]["weight"] == 1:
            q = Quadratic(name=f"edge_{i}_{j}", weight=1, rhs=1, direction=">=")
            q.addTerms(
                coeffs=[1, 1],
                vars=[x[i], x[i]],
                values=[0, 1],
                vars2=[x[j], x[j]],
                values2=[1, 0],
            )
            model.addConstraint(q)
        elif G[i][j]["weight"] == -1:
            q = Quadratic(name=f"edge_{i}_{j}", weight=1, rhs=0, direction="<=")
            q.addTerms(
                coeffs=[1, 1],
                vars=[x[i], x[i]],
                values=[0, 1],
                vars2=[x[j], x[j]],
                values2=[1, 0],
            )
            model.addConstraint(q)
        else:
            print("no such edge! ")
            break

    model.Params.TimeLimit = 100
    model.Params.RandomSeed = 123
    sol, violated = model.optimize()
    val = 0
    for (i, j) in G.edges():
        if x[i].value != x[j].value:
            if G[i][j]["weight"] == 1:
                val += 1
            elif G[i][j]["weight"] == -1:
                val -= 1
    print(iter, n, m, val)
 ================ Now solving the problem ================ 

1 800 19176 11624

SCOPとGurobi (ver. 10) ともに100秒で打ち切ったときの実験結果を以下に示す.

folder = "../data/maxcut/"
df = pd.read_csv(folder + "maxcut_bestvalues.csv", index_col=0)
df["gap"] = df["Best Known"] / df["scop(100s)"] * 100
df
#df["gurobi"] = gurobi
#df.to_csv(folder + "maxcut_bestvalues.csv")
Name Best Known Upper bound scop(100s) n m gap gurobi
0 G1 11624 12078 11624 800 19176.0 100.000000 11624.0
1 G2 11620 12084 11620 800 19176.0 100.000000 11602.0
2 G3 11622 12077 11620 800 19176.0 100.017212 11622.0
3 G4 11646 * 11646 800 19176.0 100.000000 11641.0
4 G5 11631 * 11630 800 19176.0 100.008598 11631.0
5 G6 2178 * 2172 800 19176.0 100.276243 2178.0
6 G7 2003 * 1999 800 19176.0 100.200100 1999.0
7 G8 2003 * 2005 800 19176.0 99.900249 2005.0
8 G9 2048 * 2051 800 19176.0 99.853730 2054.0
9 G10 1994 * 2000 800 19176.0 99.700000 1992.0
10 G11 564 627 564 800 19176.0 100.000000 564.0
11 G12 556 621 556 800 19176.0 100.000000 556.0
12 G13 580 645 578 800 19176.0 100.346021 582.0
13 G14 3060 3187 3049 800 4667.0 100.360774 3037.0
14 G15 3049 3169 3034 800 4667.0 100.494397 3020.0
15 G16 3045 3172 3036 800 4667.0 100.296443 3032.0
16 G17 3043 * 3030 800 4667.0 100.429043 3020.0
17 G18 988 * 984 800 4667.0 100.406504 943.0
18 G19 903 * 898 800 4667.0 100.556793 883.0
19 G20 941 * 936 800 4667.0 100.534188 882.0
20 G21 931 * 929 800 4667.0 100.215285 901.0
21 G22 13346 14123 13304 2000 19990.0 100.315695 13224.0
22 G23 13317 14129 13314 2000 19990.0 100.022533 13200.0
23 G24 13314 14131 13306 2000 19990.0 100.060123 13160.0
24 G25 13326 * 13294 2000 19990.0 100.240710 13115.0
25 G26 13314 * 13286 2000 19990.0 100.210748 13155.0
26 G27 3318 * 3296 2000 19990.0 100.667476 3183.0
27 G28 3285 * 3245 2000 19990.0 101.232666 3113.0
28 G29 3389 * 3356 2000 19990.0 100.983313 3230.0
29 G30 3403 * 3374 2000 19990.0 100.859514 3249.0
30 G31 3288 * 3250 2000 19990.0 101.169231 3088.0
31 G32 1398 1560 1386 2000 4000.0 100.865801 1410.0
32 G33 1376 1537 1368 2000 4000.0 100.584795 1382.0
33 G34 1372 1541 1364 2000 4000.0 100.586510 1384.0
34 G35 7670 8000 7613 2000 NaN 100.748719 7589.0
35 G36 7660 7996 7587 2000 NaN 100.962172 7575.0
36 G37 7666 8009 7609 2000 NaN 100.749113 7554.0
37 G38 7681 * 7618 2000 NaN 100.826989 7581.0
38 G39 2395 * 2358 2000 NaN 101.569126 2289.0
39 G40 2387 * 2341 2000 NaN 101.964972 2257.0
40 G41 2398 * 2331 2000 NaN 102.874303 2259.0
41 G42 2469 * 2422 2000 NaN 101.940545 2338.0
42 G43 6659 7027 6656 2000 NaN 100.045072 6567.0
43 G44 6648 7022 6649 1000 NaN 99.984960 6618.0
44 G45 6652 7020 6652 1000 NaN 100.000000 6612.0
45 G46 6645 * 6641 1000 NaN 100.060232 6606.0
46 G47 6656 * 6648 1000 NaN 100.120337 6608.0
47 G48 6000 * 6000 1000 NaN 100.000000 6000.0
48 G49 6000 * 6000 1000 NaN 100.000000 6000.0
49 G50 5880 * 5880 1000 NaN 100.000000 5880.0
50 G51 3846 * 3823 1000 NaN 100.601622 3808.0
51 G52 3849 * 3823 1000 NaN 100.680094 3800.0
52 G53 3846 * 3825 1000 5914.0 100.549020 3814.0
53 G54 3846 * 3828 1000 5916.0 100.470219 3811.0