Search code examples
pythonlinear-programmingscip

Too few dual values in linear model


There should be a dual solution associated with each constraint.

I have len(E)+len(D) constraints which equals 84+20 = 104

But I have only 23 dual values?

My input data is below:

#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)

My model is below

# 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.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)

This is the output result:

[49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, 49864.0, -0.0, 5168.0, 9840.0]


Solution

  • I solved it adding the three lines of code below to the model:

    model.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
    model.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
    model.disablePropagation()
    

    Reason: The model is presolved by some build in heuristics and cuts before being delivered to the LP solver, which is the reason it gave the wrong output. By setting those 'parameters' to off, the model is given directly to the LP solver, and solved correctly.