99 lines
3.0 KiB
Python
99 lines
3.0 KiB
Python
'''
|
|
gap is not constant.
|
|
|
|
OPT(IP): 3.0
|
|
OPT(LP): 0.6000000000000003
|
|
gap: 4.999999999999997
|
|
w: {(0, 2): 7, (0, 3): 6, (1, 2): 3, (1, 3): 9}
|
|
c: {(0, 2): 8, (0, 3): 9, (1, 2): 10, (1, 3): 3}
|
|
b: 8.0
|
|
edges: [(0, 2), (0, 3), (1, 2), (1, 3)]
|
|
IP solution: [('x[0,2]', 0.0), ('x[0,3]', 0.0), ('x[1,2]', 1.0), ('x[1,3]', 0.0), ('y[0,2]', 0.0), ('y[0,3]', 0.0), ('y[1,2]', 0.0), ('y[1,3]', 1.0)]
|
|
LP solution: [('x[0,2]', 0.0), ('x[0,3]', 0.0), ('x[1,2]', 0.20000000000000012), ('x[1,3]', 0.0), ('y[0,2]', 0.39999999999999997), ('y[0,3]', 0.39999999999999997), ('y[1,2]', 0.0), ('y[1,3]', 0.39999999999999997)]
|
|
'''
|
|
|
|
from random import randint
|
|
|
|
from sage.all import *
|
|
from sage.matroids.all import *
|
|
|
|
import gurobipy as gp
|
|
from gurobipy import GRB
|
|
|
|
def interdiction(G, w, c, b, vartype=GRB.BINARY):
|
|
"""
|
|
Interdiction problem on a graph G with cut weights w, interdiction costs c
|
|
and budget b.
|
|
"""
|
|
model = gp.Model("mip1")
|
|
# model.Params.LogToConsole = 0
|
|
edges=G.edges(labels=False)
|
|
# Variables
|
|
x = model.addVars(edges, vtype=vartype, name="x", lb=0)
|
|
y = model.addVars(edges, vtype=vartype, name="y", lb=0)
|
|
|
|
# Objective
|
|
model.setObjective(gp.quicksum(w[e] * x[e] for e in edges),
|
|
GRB.MINIMIZE)
|
|
|
|
# Constraints
|
|
# knapsack
|
|
model.addConstr(gp.quicksum(c[e] * y[e] for e in edges) <= b)
|
|
# cut
|
|
M=Matroid(G)
|
|
for b in M.bases():
|
|
model.addConstr(gp.quicksum(x[e]+y[e] for e in b) >= 1)
|
|
|
|
# Solve
|
|
model.optimize()
|
|
solutions = []
|
|
for v in model.getVars():
|
|
solutions.append((v.varName, v.x))
|
|
|
|
return model.ObjVal,solutions
|
|
|
|
def mincut(G, c, vartype=GRB.BINARY):
|
|
model = gp.Model("mip1")
|
|
# model.Params.LogToConsole = 0
|
|
edges=G.edges(labels=False)
|
|
# Variables
|
|
x = model.addVars(edges, vtype=vartype, name="x", lb=0)
|
|
# Objective
|
|
model.setObjective(gp.quicksum(c[e] * x[e] for e in edges),
|
|
GRB.MINIMIZE)
|
|
# Constraints
|
|
M=Matroid(G)
|
|
for b in M.bases():
|
|
model.addConstr(gp.quicksum(x[e] for e in b) >= 1)
|
|
# Solve
|
|
model.optimize()
|
|
return model.ObjVal
|
|
|
|
## enumerate all graphs
|
|
lb = 4
|
|
ub = 6
|
|
|
|
f = lambda G: G.is_connected()
|
|
for n in range(lb, ub + 1):
|
|
for G in filter(f, graphs(n)):
|
|
for _ in range(10): # use 10 random weights
|
|
# generate w,c,b
|
|
edges=G.edges(labels=False)
|
|
w = {e: randint(1, 10) for e in edges}
|
|
c = {e: randint(1, 10) for e in edges}
|
|
b = mincut(G, c)-randint(1, 5)
|
|
if b < 0: continue
|
|
|
|
IP,solIP=interdiction(G, w, c, b, GRB.BINARY)
|
|
LP,solLP=interdiction(G, w, c, b, GRB.CONTINUOUS)
|
|
print(f"OPT(IP): {IP}")
|
|
print(f"OPT(LP): {LP}")
|
|
print(f"gap: {IP/LP}")
|
|
if IP/LP > 4:
|
|
print(f"w: {w}")
|
|
print(f"c: {c}")
|
|
print(f"b: {b}")
|
|
print("edges:", G.edges(labels=False))
|
|
print("IP solution:", solIP)
|
|
print("LP solution:", solLP)
|
|
exit() |