Я пытаюсь реализовать ветвление ограничения (Райан и Фостер, 1981) в Python, используя SCIP
. Для этого мне нужно получить значения переменных в каждом ограничении.
Например, в ограничении "model.addCons (quicksum (a [e] [r] [d] * x [e, r] для e в диапазоне (len (E)) для r в диапазоне (len (R)) )> = D [d]) "Я хотел бы получить матрицу x (x [e, r]) для данного d.
Мой код указан ниже:
#Input data
from pyscipopt import Model, quicksum
def parameter_a(E,R,D):
a = [[[0 for d in range(len(D))] for r in range(len(R))] for e in range(len(E))]
for e in range(len(E)):
for d in range(len(D)):
for r in range(len(R)):
for idx_r in range(len(R[r])):
if R[r][idx_r] == 1:
if d == idx_r:
a[e][r][d] = 1
return a
def cost_matrix(E,R, arcs):
c = [[0 for r in range(len(R))] for e in range(len(E))]
for r in range(len(R)):
for e in range(len(E)):
c[e][r] = Path_Cost(E[e],R[r], arcs)
return c
E = ['nurse1', 'nurse2', 'nurse3', 'nurse4', 'nurse5', 'nurse6', 'nurse7', 'nurse8', 'nurse9', 'nurse10', 'nurse11', 'nurse12', 'nurse13', 'nurse14', 'nurse15', 'nurse16', 'nurse17', 'nurse18', 'nurse19', 'nurse20']
D = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
R = [[1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0], [0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0], [0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1]]
a = parameter_a(E,R,D)
arcs = [(0, 3, 1712.0), (0, 4, 1816.0), (0, 5, 1888.0), (1, 4, 1920.0), (1, 5, 1992.0), (2, 5, 2064.0), (3, 6, 1712.0), (3, 7, 1816.0), (3, 8, 1888.0), (4, 7, 1920.0), (4, 8, 1992.0), (5, 8, 2064.0), (6, 9, 1712.0), (6, 10, 1816.0), (6, 11, 1888.0), (7, 10, 1920.0), (7, 11, 1992.0), (8, 11, 2064.0), (9, 12, 1712.0), (9, 13, 1816.0), (9, 14, 1888.0), (10, 13, 1920.0), (10, 14, 1992.0), (11, 14, 2064.0), (12, 15, 1712.0), (12, 16, 2044.0), (12, 17, 2172.0), (13, 16, 1920.0), (13, 17, 2276.0), (14, 17, 2064.0), (15, 18, 2232.0), (15, 19, 2304.0), (15, 20, 2432.0), (16, 19, 2376.0), (16, 20, 2504.0), (17, 20, 2632.0), (18, 21, 2232.0), (18, 22, 2076.0), (18, 23, 2148.0), (19, 22, 2376.0), (19, 23, 2220.0), (20, 23, 2632.0), (21, 24, 1712.0), (21, 25, 1816.0), (21, 26, 1888.0), (22, 25, 1920.0), (22, 26, 1992.0), (23, 26, 2064.0), (24, 27, 1712.0), (24, 28, 1816.0), (24, 29, 1888.0), (25, 28, 1920.0), (25, 29, 1992.0), (26, 29, 2064.0), (27, 30, 1712.0), (27, 31, 1816.0), (27, 32, 1888.0), (28, 31, 1920.0), (28, 32, 1992.0), (29, 32, 2064.0), (30, 33, 1712.0), (30, 34, 1816.0), (30, 35, 1888.0), (31, 34, 1920.0), (31, 35, 1992.0), (32, 35, 2064.0), (33, 36, 1712.0), (33, 37, 2044.0), (33, 38, 2172.0), (34, 37, 1920.0), (34, 38, 2276.0), (35, 38, 2064.0), (36, 39, 2232.0), (36, 40, 2304.0), (36, 41, 2432.0), (37, 40, 2376.0), (37, 41, 2504.0), (38, 41, 2632.0), (39, 42, 2232.0), (39, 43, 2076.0), (39, 44, 2148.0), (40, 43, 2376.0), (40, 44, 2220.0), (41, 44, 2632.0), (42, 45, 1712.0), (42, 46, 1816.0), (42, 47, 1888.0), (43, 46, 1920.0), (43, 47, 1992.0), (44, 47, 2064.0), (45, 48, 1712.0), (45, 49, 1816.0), (45, 50, 1888.0), (46, 49, 1920.0), (46, 50, 1992.0), (47, 50, 2064.0), (48, 51, 1712.0), (48, 52, 1816.0), (48, 53, 1888.0), (49, 52, 1920.0), (49, 53, 1992.0), (50, 53, 2064.0), (51, 54, 1712.0), (51, 55, 1816.0), (51, 56, 1888.0), (52, 55, 1920.0), (52, 56, 1992.0), (53, 56, 2064.0), (54, 57, 1712.0), (54, 58, 2044.0), (54, 59, 2172.0), (55, 58, 1920.0), (55, 59, 2276.0), (56, 59, 2064.0), (57, 60, 2232.0), (57, 61, 2304.0), (57, 62, 2432.0), (58, 61, 2376.0), (58, 62, 2504.0), (59, 62, 2632.0), (60, 63, 2232.0), (60, 64, 2076.0), (60, 65, 2148.0), (61, 64, 2376.0), (61, 65, 2220.0), (62, 65, 2632.0), (63, 66, 1712.0), (63, 67, 1816.0), (63, 68, 1888.0), (64, 67, 1920.0), (64, 68, 1992.0), (65, 68, 2064.0), (66, 69, 1712.0), (66, 70, 1816.0), (66, 71, 1888.0), (67, 70, 1920.0), (67, 71, 1992.0), (68, 71, 2064.0), (69, 72, 1712.0), (69, 73, 1816.0), (69, 74, 1888.0), (70, 73, 1920.0), (70, 74, 1992.0), (71, 74, 2064.0), (72, 75, 1712.0), (72, 76, 1816.0), (72, 77, 1888.0), (73, 76, 1920.0), (73, 77, 1992.0), (74, 77, 2064.0), (75, 78, 1712.0), (75, 79, 2044.0), (75, 80, 2172.0), (76, 79, 1920.0), (76, 80, 2276.0), (77, 80, 2064.0), (78, 81, 2232.0), (78, 82, 2304.0), (78, 83, 2432.0), (79, 82, 2376.0), (79, 83, 2504.0), (80, 83, 2632.0)]
c = cost_matrix(E,R, arcs)
# NRP Model
def NRP():
model = Model("NRP")
# Variables
x = {}
for e in range(len(E)):
for r in range(len(R)):
x[e,r] = model.addVar(vtype="C", name='x'+str(e)+'.'+str(r))
# Constraints
for e in range(len(E)):
model.addCons(quicksum(x[e,r] for r in range(len(R))) == 1)
for d in range(len(D)):
model.addCons(quicksum(a[e][r][d]*x[e,r] for e in range(len(E)) for r in range(len(R))) >= D[d])
# Objective
objective = quicksum(c[e][r]*x[e,r] for e in range(len(E)) for r in range(len(R)))
model.setObjective(objective, "minimize")
return model
def Solve(model, print_satement=False):
model.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
model.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
model.disablePropagation()
model.optimize()
optimal = 0
primal = []
dual = []
if model.getStatus() == "optimal":
optimal = model.getObjVal()
if print_satement == True:
print("Optimal value:", model.getObjVal())
for v in model.getVars():
primal += [model.getVal(v)]
if print_satement == True:
print(v.name, " = ", model.getVal(v))
for c in model.getConss():
dual += [model.getDualsolLinear(c)]
if print_satement == True:
print(c.name, model.getSlack(c), model.getDualsolLinear(c))
elif print_satement == True:
print("Problem could not be solved to optimality")
return model, optimal, primal, dual
model = NRP()
model, optimal, primal, dual = Solve(model, print_satement = False)
print(dual)