In [1]:
from gurobipy import *
In [3]:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 24 22:05:03 2018
@author: yoichi_izunaga
"""
#---packages---#
import numpy as np
from gurobipy import *
#--------------#
"""
#---read data---#
A_org: adjacency matrix of the original graph
A: agjacency matrix of the graph G with the source node(0)
n: the number of nodes of G
m: the number of directed-edges of G
"""
A = np.array(
[[0.0, 0.0, 1.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
]
)
n = len(A[0])
m = int(np.sum(A))
extA = np.c_[np.zeros((len(A)+1,1)), np.r_[ np.ones((1,len(A))), A] ]
print("n = {0:}".format(n))
print("m = {0:}".format(m))
"""
#---construct each set---#
V: the set of nodes of G
E: the set of directed-edges of G
extV: the set of nodes of the extended graph
extE: the set of edges of the extended graph
delta: the set of tail nodes of each node
rho: the set of head nodes of each node
"""
V = list(range(1,n+1))
extV = list(range(n+1))
extE = []
for i in extV:
for j in extV:
if extA[i,j]==1.0:
extE.append((i,j))
E = []
E = [(i,j) for (i,j) in extE if i!=0 and j!=0]
delta = []
rho = []
for i in extV:
#delta.append(A[i].tolist())
delta.append(np.where(extA[i]==1.0)[0].tolist())
rho.append(np.where(extA[:,i]==1.0)[0].tolist())
print("元の頂点集合 = {0:}".format(V))
print("頂点集合 = {0:}".format(extV))
print("元の枝集合 = {0:}".format(E))
print("枝集合 = {0:}".format(extE))
print("tail = {0:}".format(delta))
print("head = {0:}".format(rho))
AP = Model("TreeGen")
y = {}
z = {}
l = {}
for i in extV:
y[i] = AP.addVar(vtype="B", name="y(%s)"%i)
l[i] = AP.addVar(vtype="C", name="l(%s)"%i)
for (i,j) in extE:
z[i,j] = AP.addVar(vtype="B", name="z(%s,%s)"%(i,j))
AP.update()
#ソースノードから1本だけ出る
AP.addConstr( quicksum(z[0,j] for j in delta[0])==1 )
#とりあえず頂点数を固定
AP.addConstr( quicksum(y[i] for i in V)==5 )
#頂点数-1=枝数
AP.addConstr( quicksum(y[i] for i in V)-1 == quicksum(z[i,j] for (i,j) in E) )
#Tree上の頂点の入次数は1
for i in V:
AP.addConstr( y[i] == quicksum(z[j,i] for j in rho[i]) )
# tourを禁止する制約
for (i,j) in E:
AP.addConstr( l[i]+1 <= l[j]+n*(1-z[i,j]) )
AP.setObjective(quicksum(z[i,j] for (i,j) in E), GRB.MINIMIZE)
AP.update()
AP.write("test.lp")
AP.optimize()
sol_y = [ y[i].X for i in V ]
sol_z = [ z[i,j].X for (i,j) in E ]
#sol_l = [ l[i].X for i in range(n) ]
T = [ (i,j) for (i,j) in E if z[i,j].X==1]
print("y = \n", sol_y)
print("z = \n", sol_z)
#print("l = \n", sol_l)
print(T)
In [55]:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 24 22:05:03 2018
@author: yoichi_izunaga
"""
#---packages---#
import numpy as np
from gurobipy import *
#--------------#
"""
#---read data---#
A_org: adjacency matrix of the original graph
A: agjacency matrix of the graph G with the source node(0)
n: the number of nodes of G
m: the number of directed-edges of G
"""
D = np.array([[0,1,1,0,0,1],
[0,0,1,0,0,0],
[0,1,0,0,0,0],
[1,0,0,0,0,1],
[1,0,1,0,0,0],
[0,1,1,0,1,0]])
n = len(A[0])
m = int(np.sum(A))
print("n = {0:}".format(n))
print("m = {0:}".format(m))
AP = Model("rankability")
x = {}
y = {}
for i in range(n):
for j in range(n):
x[i,j] = AP.addVar(vtype="C",ub=1,name="x(%s,%s)"%(i,j))
y[i,j] = AP.addVar(vtype="C",ub=1,name="y(%s,%s)"%(i,j))
AP.update()
In [59]:
for i in range(n):
for j in range(n):
AP.addConstr(x[i,j] + y[i,j] <= 1)
AP.addConstr(x[i,j] <= 1 - D[i,j])
AP.addConstr(y[i,j] <= D[i,j])
for j in range(1,n):
for i in range(0,j):
AP.addConstr(D[i,j] + x[i,j] - y[i,j] + D[j,i] + x[j,i] - y[j,i] == 1)
for i in range(n):
for j in range(n):
for k in range(n):
if j != i and k != j and k != i:
AP.addConstr(D[i,j] + x[i,j] - y[i,j] + D[j,k] + x[j,k] - y[j,k] + D[k,i] + x[k,i] - y[k,i] <= 2)
smartk = n*n
AP.addConstr(quicksum(x[i,j] for i in range(n) for j in range(n)) + quicksum(y[i,j] for i in range(n) for j in range(n)) <= smartk)
AP.update()
AP.setObjective(quicksum(x[i,j] for i in range(n) for j in range(n)) + quicksum(y[i,j] for i in range(n) for j in range(n)),GRB.MINIMIZE)
AP.update()
AP.write('rankability.lp')
#LP = AP.relax()
AP.Params.Method = 2
AP.Params.Crossover = 0
AP.update()
AP.optimize()
print(AP.Status)
In [61]:
sol_x = [round(100*x[i,j].X)/100 for i in range(n) for j in range(n)]
sol_y = [round(100*y[i,j].X)/100 for i in range(n) for j in range(n)]
sol_x = np.reshape(sol_x,(n,n))
sol_y = np.reshape(sol_y,(n,n))
print(sol_x)
print(sol_y)
In [54]:
dir(x[0,0].Xn)
In [ ]: