Movatterモバイル変換

[0]ホーム
URL:

グラフ分割問題

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

準備

ここでは，付録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 graphParameters:

- 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 problemParameters:

- 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.__dataL=[iforiinVifx[i].X>=0.5]R=[iforiinVifx[i].X<0.5]
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)Thread count: 8 physical cores, 16 logical processors, using up to 16 threadsOptimize a model with 55 rows, 37 columns and 172 nonzerosModel fingerprint: 0x44488fa2Variable types: 0 continuous, 37 integer (37 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 13.0000000Presolve time: 0.00sPresolved: 55 rows, 37 columns, 172 nonzerosVariable types: 0 continuous, 37 integer (37 binary)Root relaxation: objective 0.000000e+00, 16 iterations, 0.00 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time     0     0    0.00000    0   10   13.00000    0.00000   100%     -    0sH    0     0                      12.0000000    0.00000   100%     -    0sH    0     0                       9.0000000    4.66667  48.1%     -    0s     0     0    8.00000    0   31    9.00000    8.00000  11.1%     -    0sCutting planes:  Gomory: 5  MIR: 2  StrongCG: 1  Zero half: 2  Mod-K: 2  RLT: 19Explored 1 nodes (78 simplex iterations) in 0.02 secondsThread count was 16 (of 16 available processors)Solution count 3: 9 12 13 Optimal solution found (tolerance 1.00e-04)Best objective 9.000000000000e+00, best bound 9.000000000000e+00, gap 0.0000%Opt.value= 9.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)inEif(iinLandjinR)or(iinRandjinL)]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 toEquipart: 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.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)

タブーサーチ

タブーサーチ(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 goingfrom 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[source]

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

Execute a tabu search run.

rndseed=2random.seed(rndseed)n=100pos={i:(random.random(),random.random())foriinrange(n)}G=nx.random_geometric_graph(n,0.2,pos=pos)nodes=G.nodes()adj=[set([])foriinnodes]for(i,j)inG.edges():adj[i].add(j)adj[j].add(i)max_iterations=1000tabulen=len(nodes)/2sol=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.0final 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=[iforiinnodesifsol[i]==1]R=[iforiinnodesifsol[i]==0]edgelist=[(i,j)for(i,j)inG.edges()if(iinLandjinR)or(iinRandjinL)]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 estimatedtemperature 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 probabilitysizefactor=16# for det. # tentatives at each temp.tempfactor=0.95# cooling ratiofreezelim=5# max number of iterations with less that minpercent acceptancesminpercent=0.02# fraction of accepted moves for being not frozenalpha=3.0# imballance factorN=len(nodes)# neighborhood sizeL=N*sizefactor# number of tentatives at current temperatureprint("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 solutionz,bal,s,d=evaluate(nodes,adj,sol,alpha)print("initial partition: z =",z)print(sol)print()LOG=Falsesol,z=annealing(nodes,adj,sol,initprob,L,tempfactor,freezelim,minpercent,alpha,report=False)zp,bal,sp,dp=evaluate(nodes,adj,sol,alpha)assertabs(zp-z)<1.0e-9# floating point approx.equalityprint()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.0initial partition: z = 226.0[0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0]final solution: z= 33.0[1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1]
L=[iforiinnodesifsol[i]==1]R=[iforiinnodesifsol[i]==0]edgelist=[(i,j)for(i,j)inG.edges()if(iinLandjinR)or(iinRandjinL)]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)
#     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=[iforiinnodesifsol[i]==1]R=[iforiinnodesifsol[i]==0]edgelist=[(i,j)for(i,j)inG.edges()if(iinLandjinR)or(iinRandjinL)]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)とよばれ，様々な制約が付加された問題が考えられている．

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

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

METISで最適化する関数 metis

引数

  • $G$: networkXの無向グラフのインスタンス
  • npartition: 分割数を表す2以上の整数値（既定値は2）
  • ufactor: 分割に含まれる点数の非均衡度の上限を与えるパラメータ（既定値は30であり，30\%の非均衡を上限とする．）

metis[source]

metis(G:Graph,npartition:int=2,ufactor:int=30)

metis関数の使用例

rndseed=2random.seed(rndseed)n=10000pos={i:(random.random(),random.random())foriinrange(n)}G=nx.random_geometric_graph(n,0.05,pos=pos)#枝に重みを付与foru,vinG.edges():G[u][v]["weight"]=random.randint(1,1)#点に重みを付与foruinG.nodes():G.nodes[u]["weight"]=random.randint(1,100)nx.draw(G,pos=pos,node_size=0.3)
npartition=5#分割数ufactor=30#不均衡パラメータsol=metis(G,npartition,ufactor)num_nodes=np.zeros(npartition,int)foriinsol:num_nodes[sol[i]]+=G.nodes[i]["weight"]color_list=[sol[i]foriinG.nodes()]nx.draw_networkx_nodes(G,pos,node_color=color_list,node_size=5)# 枝の描画edge_colors=[]edge_widths=[]obj=0foru,vinG.edges():ifsol[u]!=sol[v]:# 両端ノードの色が異なる場合obj+=G[u][v]["weight"]edge_colors.append("r")# 赤色で描画edge_widths.append(1)# 太さを2に設定else:edge_colors.append("k")# 黒色で描画edge_widths.append(0.1)# 太さを1に設定nx.draw_networkx_edges(G,pos,edge_color=edge_colors,width=edge_widths)print(obj,num_nodes)
16899 [103106  99953 101341  98676 102434]

最大カット問題

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

無向グラフ $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}$$
defmaxcut(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={}foriinV:x[i]=model.addVar(vtype="B",name=f"x({i})")for(i,j)inE:y[i,j]=model.addVar(vtype="B",name=f"y({i},{j})")model.update()for(i,j)inE: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)inE),GRB.MAXIMIZE)model.update()model.__data=xreturnmodel

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}$$
defmaxcut_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={},{},{}foriinV:x[i]=model.addVar(vtype="B",name=f"x({i})")for(i,j)inE: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)inE: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)inE),GRB.MAXIMIZE)#model.update()#model.__data = xreturnmodeldefmake_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)foriinVforjinVifi<jandrandom.random()<prob]returnV,E
#GRB = MDOfromgurobipyimportModel,quicksum,GRBrandom.seed(3)V,E=make_data(100,0.1)model=maxcut(V,E)model.optimize()status=model.Statusifstatus==GRB.Status.OPTIMAL:print("Opt.value=",model.ObjVal)obj_lin=model.ObjValx=model.__dataprint([iforiinVifx[i].X>=0.5])print([iforiinVifx[i].X<0.5])else:pass
#GRB = MDOfromgurobipyimportModel,quicksum,GRBmodel=maxcut_soco(V,E)model.optimize()status=model.Statusifstatus==GRB.Status.OPTIMAL:print("Opt.value=",model.ObjVal)x=model.__dataprint([iforiinVifx[i].X>=0.5])print([iforiinVifx[i].X<0.5])# assert model.ObjVal == obj_linelse:pass

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/

fromgurobipyimportModel,quicksum,GRBfolder="../data/maxcut/"gurobi=[]foriterinrange(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=[]forrowindata[1:]:i,j,w=list(map(int,row.split()))edges.append((i,j,w))G.add_weighted_edges_from(edges)model=Model()x={}foriinG.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)inG.edges()),GRB.MAXIMIZE)model.Params.TimeLimit=10model.Params.OutputFlag=1model.optimize()print(iter,model.ObjVal)gurobi.append(model.ObjVal)

制約最適化ソルバー 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モジュールからインポートする．

fromscopimport*
folder="../data/maxcut/"foriterinrange(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=[]forrowindata[1:]:i,j,w=list(map(int,row.split()))edges.append((i,j,w))G.add_weighted_edges_from(edges)model=Model()x={}foriinG.nodes():x[i]=model.addVariable(name=f"x[{i}]",domain=[0,1])for(i,j)inG.edges():ifG[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)elifG[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! ")breakmodel.Params.TimeLimit=100model.Params.RandomSeed=123sol,violated=model.optimize()val=0for(i,j)inG.edges():ifx[i].value!=x[j].value:ifG[i][j]["weight"]==1:val+=1elifG[i][j]["weight"]==-1:val-=1print(iter,n,m,val)

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)"]*100df#df["gurobi"] = gurobi#df.to_csv(folder + "maxcut_bestvalues.csv")
[8]ページ先頭
©2009-2025 Movatter.jp