Movatterモバイル変換

[0]ホーム
URL:

k-センター問題

k-センター問題に対する定式化と解法

準備

importrandomimportmathfromgurobipyimportModel,quicksum,GRB# from mypulp import Model, quicksum, GRBimportnetworkxasnximportmatplotlib.pyplotasplt

センター問題

センター問題（center problem）とは，顧客から最も近い施設への距離の「最大値」を最小にするように，グラフ内の点または枝上，または空間内の任意の点から $k$ 個の施設を選択する問題の総称であり，施設数 $k$ を指定した問題を$k$-センター問題($k$-center problem) とよぶ．

ここでは，点上に施設を配置する場合を考える．

標準定式化

$k$-メディアン問題と同様に，顧客 $i$ から施設 $j$ への距離を $c_{ij}$ とし，以下の変数を導入する．

$$ x_{ij}= \left\{ \begin{array}{ll} 1 & 顧客 i の需要が施設 j によって満たされるとき \\ 0 & それ以外のとき \end{array} \right.$$$$ y_j = \left\{ \begin{array}{ll} 1 & 施設 j を開設するとき \\ 0 & それ以外のとき \end{array} \right.$$

最も遠い施設でサービスを受ける顧客の移動費用を表す連続変数 $z$ を導入する．

上の記号および変数を用いると，$k$-センター問題は以下のように定式化できる．

$$ \begin{array}{l l l} minimize & z & \\s.t. & \sum_{j \in J} x_{ij} =1 & \forall i \in I \\ & \sum_{j \in J} c_{ij} x_{ij} \leq z & \forall i \in I \\ & \sum_{j \in J} y_{j} = k & \\ & x_{ij} \leq y_j & \forall i \in I, j \in J \\ & x_{ij} \in \{ 0,1 \} & \forall i \in I, j \in J \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array}$$
defdistance(x1,y1,x2,y2):returnmath.sqrt((x2-x1)**2+(y2-y1)**2)defmake_data(n):I=range(n)J=range(n)x=[random.random()foriinrange(n)]y=[random.random()foriinrange(n)]c={}foriinI:forjinJ:c[i,j]=distance(x[i],y[i],x[j],y[j])returnI,J,c,x,ydefkcenter(I,J,c,k):model=Model("k-center")z=model.addVar(vtype="C")x,y={},{}forjinJ:y[j]=model.addVar(vtype="B")foriinI:x[i,j]=model.addVar(vtype="B")model.update()foriinI:model.addConstr(quicksum(x[i,j]forjinJ)==1)model.addConstr(quicksum(c[i,j]*x[i,j]forjinJ)<=z)forjinJ:model.addConstr(x[i,j]<=y[j])model.addConstr(quicksum(y[j]forjinJ)==k)model.setObjective(z,GRB.MINIMIZE)model.update()model.__data=x,yreturnmodel
random.seed(67)n=30I,J,c,x_pos,y_pos=make_data(n)k=5model=kcenter(I,J,c,k)model.optimize()x,y=model.__datafacilities=[jforjinyify[j].X>0.5]edges=[(i,j)for(i,j)inxifx[i,j].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 961 rows, 931 columns and 3630 nonzerosModel fingerprint: 0x8230ce33Variable types: 1 continuous, 930 integer (930 binary)Coefficient statistics:  Matrix range     [3e-02, 1e+00]  Objective range  [1e+00, 1e+00]  Bounds range     [1e+00, 1e+00]  RHS range        [1e+00, 5e+00]Presolve time: 0.01sPresolved: 961 rows, 931 columns, 3630 nonzerosVariable types: 1 continuous, 930 integer (930 binary)Found heuristic solution: objective 0.7069622Found heuristic solution: objective 0.6501211Found heuristic solution: objective 0.6332376Root relaxation: objective 1.998137e-01, 507 iterations, 0.01 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time     0     0    0.19981    0  228    0.63324    0.19981  68.4%     -    0sH    0     0                       0.5951071    0.19981  66.4%     -    0sH    0     0                       0.4970149    0.19981  59.8%     -    0s     0     0    0.22202    0  222    0.49701    0.22202  55.3%     -    0s     0     0    0.22202    0  228    0.49701    0.22202  55.3%     -    0sH    0     0                       0.4358186    0.22202  49.1%     -    0sH    0     0                       0.4298363    0.22202  48.3%     -    0sH    0     0                       0.4044511    0.22221  45.1%     -    0s     0     0    0.22952    0  193    0.40445    0.22952  43.3%     -    0sH    0     0                       0.3950400    0.22952  41.9%     -    0s     0     0    0.23666    0  169    0.39504    0.23666  40.1%     -    0s     0     0    0.24335    0  161    0.39504    0.24335  38.4%     -    0s     0     0    0.24335    0  129    0.39504    0.24335  38.4%     -    0sH    0     0                       0.3700418    0.24335  34.2%     -    0sH    0     0                       0.3455534    0.24335  29.6%     -    0s     0     0    0.26065    0  120    0.34555    0.26065  24.6%     -    0s     0     0    0.26445    0  112    0.34555    0.26445  23.5%     -    0s     0     0    0.27611    0  119    0.34555    0.27611  20.1%     -    0s     0     0    0.29390    0  101    0.34555    0.29390  14.9%     -    0s     0     0    0.29594    0  119    0.34555    0.29594  14.4%     -    0s     0     0    0.29740    0  124    0.34555    0.29740  13.9%     -    0s     0     0    0.29740    0  124    0.34555    0.29740  13.9%     -    0sH    0     0                       0.3352990    0.29740  11.3%     -    0s     0     0    0.30794    0   32    0.33530    0.30794  8.16%     -    0sH    0     0                       0.3336685    0.30794  7.71%     -    0s     0     0     cutoff    0         0.33367    0.33367  0.00%     -    0sCutting planes:  Gomory: 6  MIR: 5  Zero half: 1  Mod-K: 1  RLT: 6Explored 1 nodes (3018 simplex iterations) in 0.27 secondsThread count was 16 (of 16 available processors)Solution count 10: 0.333668 0.335299 0.345553 ... 0.595107Optimal solution found (tolerance 1.00e-04)Best objective 3.336684909605e-01, best bound 3.336684909605e-01, gap 0.0000%
defdraw_center(facilities,edges,x_pos,y_pos):G=nx.Graph()facilities=set(facilities)unused=set(jforjinJifjnotinfacilities)client=set(iforiinIifinotinfacilitiesandinotinunused)G.add_nodes_from(facilities)G.add_nodes_from(client)G.add_nodes_from(unused)for(i,j)inedges:G.add_edge(i,j)position={}foriinrange(len(x_pos)):position[i]=(x_pos[i],y_pos[i])nx.draw(G,position,with_labels=False,node_color="b",nodelist=facilities)nx.draw(G,position,with_labels=False,node_color="c",nodelist=unused,node_size=50)nx.draw(G,position,with_labels=False,node_color="g",nodelist=client,node_size=50)draw_center(facilities,edges,x_pos,y_pos)

2分探索を用いた解法

標準定式化はmin-max型の目的関数をもつので大規模問題例には不向きである．ここでは，2分探索を用いた解法を考える．

顧客から施設までの距離が $\theta$ 以下である枝だけから構成されるグラフ $G_{\theta}=(V,E_{\theta})$ を考える．点の部分集合 $S(\subseteq V)$ が与えられたとき，すべての点 $i(\in V)$ が $S$ 内の少なくとも1つの点に隣接するとき$S$ を点被覆(vertex cover)と呼ぶ．

グラフ $G_{\theta}$ 上で $|S|=k$ の被覆が存在するなら，$k$-センター問題の最適値が $\theta$ 以下であることが言える．

施設 $j$ を開設するとき $1$ となる $0$-$1$ 変数 $y_j$ は，ここでは点の部分集合 $S$ に含まれるとき $1$ となる変数と見なすことができる．さらに以下の変数を導入する．$$ z_i = \left\{ \begin{array}{ll} 1 & 点 i が S 内の点に隣接しない(被覆されない） \\ 0 & それ以外のとき \end{array} \right.$$

グラフ $G_{\theta}$ の隣接行列（点 $i,j$ が隣接しているとき $1$，それ以外のとき $0$ の要素をもつ行列）を $[a_{ij}]$ としたとき， $|S|=k$ の被覆が存在するか否かを判定する問題は，以下の被覆立地問題(covering location problem)として定式化できる．

$$\begin{array}{l l l} minimize & \sum_{i \in I} z_{i} & \\s.t. & \sum_{j \in J} a_{ij} y_{j} +z_{i} =1 & \forall i \in I \\ & \sum_{j \in J} y_{j} = k & \\ & z_{i} \in \{ 0,1 \} & \forall i \in I \\ & y_j \in \{ 0,1 \} & \forall j \in J\end{array}$$

目的関数は，被覆されない顧客の数を最小化することを表す．最初の制約は，顧客 $i$ が $S$ 内の点にグラフ $G_{\theta}$ 上で隣接するか，そうでない場合には変数 $z_i$ が $1$ になることを規定する．2番目の制約は，開設された施設の数が $k$ 個であることを規定する．

この問題を $G_{\theta}$ 上の $k$点被覆問題と呼ぶ．$k$点被覆問題を部分問題として用いて最適な $\theta$ （上の問題の最適値が $0$ になる最小の $\theta$）は2分探索で求めることができる．

defkcover(I,J,c,k):"""kcover -- minimize the number of uncovered    customers from k facilities.    Parameters:        - I: set of customers        - J: set of potential facilities        - c[i,j]: cost of servicing customer i from facility j        - k: number of facilities to be used    Returns a model, ready to be solved.    """model=Model("k-center")z,y,x={},{},{}foriinI:z[i]=model.addVar(vtype="B",name="z(%s)"%i)forjinJ:y[j]=model.addVar(vtype="B",name="y(%s)"%j)foriinI:x[i,j]=model.addVar(vtype="B",name="x(%s,%s)"%(i,j))model.update()foriinI:model.addConstr(quicksum(x[i,j]forjinJ)+z[i]==1,"Assign(%s)"%i)forjinJ:model.addConstr(x[i,j]<=y[j],"Strong(%s,%s)"%(i,j))model.addConstr(quicksum(y[j]forjinJ)==k,"k_center")model.setObjective(quicksum(z[i]foriinI)+0.001*quicksum(c[i,j]*x[i,j]for(i,j)inx),GRB.MINIMIZE,)model.update()model.__data=x,y,zreturnmodel
defsolve_kcenter(I,J,c,k,delta):"""solve_kcenter -- locate k facilities minimizing    distance of most distant customer.    Parameters:        I - set of customers        J - set of potential facilities        c[i,j] - cost of servicing customer i from facility j        k - number of facilities to be used        delta - tolerance for terminating bisection    Returns:        - list of facilities to be used        - edges linking them to customers    """model=kcover(I,J,c,k)x,y,z=model.__datafacilities,edges=[],[]LB=0UB=max(c[i,j]for(i,j)inc)whileUB-LB>delta:theta=(UB+LB)/2.0# print "\n\ncurrent theta:", thetaforjinJ:foriinI:ifc[i,j]>theta:x[i,j].UB=0else:x[i,j].UB=1.0model.update()model.Params.OutputFlag=0# silent modemodel.Params.Cutoff=0.1model.optimize()ifmodel.status==GRB.Status.OPTIMAL:UB=thetafacilities=[jforjinyify[j].X>0.5]edges=[(i,j)for(i,j)inxifx[i,j].X>0.5]else:# infeasibility > 0:LB=thetareturnfacilities,edges
random.seed(67)n=200I,J,c,x_pos,y_pos=make_data(n)k=5delta=1.0e-4facilities,edges=solve_kcenter(I,J,c,k,delta)print("Selected facilities:",facilities)print("Max distance from a facility to a customer: ",max([c[i,j]for(i,j)inedges]))
Selected facilities: [89, 93, 131, 133, 196]Max distance from a facility to a customer:  0.3061769941895321
draw_center(facilities,edges,x_pos,y_pos)
[8]ページ先頭
©2009-2025 Movatter.jp