Solving BPP using B&P in pyscipopt

[abb-omidi]

Dear community team,

First of all, thanks for your interesting SCIP python wrapper. I would like to use SCIP as a branch and price framework and am trying to do on the win10 using CMake. It works fine but, as I don't have a good background in the C/C++ programming, it is a bit hard to understand its functions. I read other issues about that and understand it's not possible for some examples like test_coloring. (this and this issues.).
I'm wondering if, is there BPP or VRP example in which implementing BP/BCP in python?

Regards
Abbas

[fserra]

Hi Abbas

Sorry, but currently there is no such example in python, only the ones available in SCIP (https://scipopt.org/doc/html/BINPACKING_MAIN.php and https://scipopt.org/doc/html/VRP_MAIN.php)
The coloring was translated from https://scipopt.org/doc/html/COLORING_MAIN.php but back then we were trying to test some extensions to the interface, thus, sadly, it doesn't work with the current pyscipopt.

I might spend some time trying to make the coloring work, but I don't promise anything.

[Steve8881]

Hi there

I am trying to implement a branch-and-price algorithm in Python to solve the bin packing problem, mainly based on the coloring example, the pricer example and the branching example in PySCIPOpt, and the bin packing example in SCIP. When Model.branchVar() is applied, everything works well (though the performance is not satisfying). However, when I try to implement the Ryan/Foster branching rule with a SameDiff constraint handler, a segmentation fault occurs (Segmentation fault: 11). After some debugging, the problem seems to be caused by the consactive() method and I have no idea how to do the further debugging. Could you guys give me some guidance to solve the problem?

Thank you.

Steve

PS: I have implemented the Model.repropagateNode() function by myself and used the Model.getObjVal() function to replace the SCIPgetLocalTransEstimate() function. Here is the code:

from pyscipopt import *
import numpy as np
import inspect

entering = lambda: print("IN ", inspect.stack()[1][3]) # debugging magic
leaving = lambda: print("OUT ", inspect.stack()[1][3]) # debugging magic

class SameDiff(Conshdlr):

    def consdataCreate(self, name, item1, item2, constype, node):
        entering()
        cons = self.model.createCons(self, name, stickingatnode = True)
        cons.data = {}
        cons.data['item1'] = item1
        cons.data['item2'] = item2
        cons.data['type'] = constype
        cons.data['npropagatedvars'] = 0
        cons.data['npropagations'] = 0
        cons.data['propagated'] = False
        cons.data['node'] = node
        leaving()
        return cons
    
    def checkVariable(self, cons, var, varid, nfixedvars):
        cutoff = False
        if var.getUbLocal() < 0.5:
            return nfixedvars, cutoff

        constype = cons.data['type']
        packings = self.model.data['packings']
        packingId = varid
        existitem1 = packings[packingId][cons.data['item1']]
        existitem2 = packings[packingId][cons.data['item2']]
        
        if (constype == 'same' and (existitem1 != existitem2)) or (constype == 'diff' and existitem1 and existitem2):
            infeasible, fixed = self.model.fixVar(var, 0.0)
            if infeasible:
                cutoff = True
            else:
                nfixedvars += 1
        return nfixedvars, cutoff

    def consdataFixVariables(self, cons):
        nfixedvars = 0
        cutoff = False
        v = cons.data['npropagatedvars']
        while (v < len(self.model.data['var'])) and (not cutoff):
            nfixedvars, cutoff = self.checkVariable(cons, self.model.data['var'][v], v, nfixedvars)
            v += 1
        if (cutoff):
            return SCIP_RESULT.CUTOFF
        elif (nfixedvars > 0):
            return SCIP_RESULT.REDUCEDDOM
        
    def consactive(self, constraint):
        entering()        
        if constraint.data['npropagatedvars'] != len(self.model.data['var']):
            constraint.data['propagated'] = False  
            self.model.repropagateNode(constraint.data['node'])
        leaving()

    def consdeactive(self, constraint):
        entering()
        constraint.data['npropagatedvars'] = len(self.model.data['var'])
        leaving()

    def consprop(self, constraints, nusefulconss, nmarkedconss, proptiming):
        entering()
        for c in constraints:
            if not c.data['propagated']:
                result = self.consdataFixVariables(c)
                c.data['npropagations'] += 1
                if result != SCIP_RESULT.CUTOFF:
                    c.data['propagated'] = True
                    c.data['npropagatedvars'] = len(self.model.data['var'])
                else:
                    break
        leaving()

class MyRyanFosterBranching(Branchrule):

    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):
        entering()
        
        lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars = self.model.getLPBranchCands()

        fractionalities = np.zeros([len(self.model.data['widths']), len(self.model.data['widths'])])

        for ii in range(len(lpcands)):

            solval = lpcandssol[ii]

            packingNum = self.model.data['varnames'].index(lpcands[ii].name)

            currentPacking = self.model.data['packings'][packingNum]

            con_id = np.where(np.array(currentPacking) == 1)[0]

            for i in con_id:
                for j in con_id:
                    if i == j:
                        continue
                    fractionalities[i, j] += solval

        temp = np.fmin(fractionalities, 1 - fractionalities)

        x_index, y_index = np.where(temp == np.max(temp))

        row1 = x_index[0]

        row2 = y_index[0]

        childsame = self.model.createChild(0.0, self.model.getObjVal())

        childdiff = self.model.createChild(0.0, self.model.getObjVal())

        conssame = self.model.data['conshdlr'].consdataCreate("some_same_name", row1, row2, "same", childsame)

        consdiff = self.model.data['conshdlr'].consdataCreate("some_diff_name", row1, row2, "diff", childdiff)

        self.model.addConsNode(childsame, conssame)
        
        self.model.addConsNode(childdiff, consdiff)

        leaving()

        return {"result": SCIP_RESULT.BRANCHED}
        
class MyVarBranching(Branchrule):

    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):

        lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars = self.model.getLPBranchCands()

        integral = lpcands[0]

        self.model.branchVar(integral)

        return {"result": SCIP_RESULT.BRANCHED}

      
class PackingPricer(Pricer):

    # The reduced cost function for the variable pricer
    def pricerredcost(self):

        # Retrieving the dual solutions
        dualSolutions = []
        for i, c in enumerate(self.data['cons']):
            dualSolutions.append(self.model.getDualsolLinear(c))

        #################################################################
        #################################################################

        # Solve the MIP directly

        # Building a MIP to solve the subproblem
        subMIP = Model("BinPacking-Sub")

        # Setting the verbosity level to 0
        subMIP.hideOutput()

        binWidthVars = []
        varNames = []
        varBaseName = "Packing"

        # Variables for the subMIP
        for i in range(len(dualSolutions)):
            varNames.append(varBaseName + "_" + str(i))
            binWidthVars.append(subMIP.addVar(varNames[i], vtype = "B", obj = -1.0 * dualSolutions[i]))

        # Adding the knapsack constraint
        subMIP.addCons(quicksum(w*v for (w,v) in zip(self.model.data['widths'], binWidthVars)) <= self.model.data['rollLength'])

        for ii in range(len(self.model.data['var'])):
            if self.model.data['var'][ii].getUbLocal() < 0.5:
                currentPacking = self.model.data['packings'][ii]
                subMIP.addCons(quicksum(w*v for (w,v) in zip(currentPacking, binWidthVars)) <= np.sum(currentPacking) - 1)


        # Solving the subMIP to generate the most negative reduced cost packing
        subMIP.optimize()

        objval = 1 + subMIP.getObjVal()

        #################################################################
        #################################################################
        

        # Adding the column to the master problem
        if objval < -1e-08:
            currentNumVar = len(self.model.data['var'])

            # Creating new var; must set pricedVar to True
            newVar = self.model.addVar("NewPacking_" + str(currentNumVar), vtype = "I", obj = 1.0, pricedVar = True)

            self.model.data['varnames'].append("NewPacking_" + str(currentNumVar))

            # Adding the new variable to the constraints of the master problem
            newPacking = []
            for i, c in enumerate(self.data['cons']):
                coeff = round(subMIP.getVal(binWidthVars[i]))
                self.model.addConsCoeff(c, newVar, coeff)

                newPacking.append(coeff)
            
            self.model.data['packings'].append(newPacking)
            self.model.data['var'].append(newVar)

        return {'result':SCIP_RESULT.SUCCESS}

    # The initialisation function for the variable pricer to retrieve the transformed constraints of the problem
    def pricerinit(self):
        for i, c in enumerate(self.data['cons']):
            self.data['cons'][i] = self.model.getTransformedCons(c)


def test_binpacking():

    s = Model()

    s.setPresolve(SCIP_PARAMSETTING.OFF)

    pricer = PackingPricer(priority = 5000000)

    s.includePricer(pricer, "binPackingPricer", "Pricer to identify new packings")

    conshdlr = SameDiff()

    s.includeConshdlr(conshdlr, "SameDiff", "Generate the SameDiff constraints",
                          chckpriority=2000000, propfreq=1)

    my_branchrule = MyRyanFosterBranching(s)

    s.includeBranchrule(my_branchrule, "MyRyanFosterBranching", "MyRyanFosterBranching",
                    priority=10000000, maxdepth= -1, maxbounddist=1)

    # item widths
    widths = [71,80,76,83,41,78,70,23,42,87,43,84,60,55,49]
    # width demand
    demand = [1] * 15 
    # roll length
    rollLength = 150
    assert len(widths) == len(demand)

    # adding the initial variables
    PackingVars = []
    varNames = []
    varBaseName = "Packing"
    packings = []

    for i in range(len(widths)):
        varNames.append(varBaseName + "_" + str(i))
        PackingVars.append(s.addVar(varNames[i], vtype = 'B', obj = 1.0))

    # adding a linear constraint for the knapsack constraint
    demandCons = []
    for i in range(len(widths)):

        numWidthsPerRoll = 1 #float(int(rollLength/widths[i]))
        demandCons.append(s.addCons(numWidthsPerRoll*PackingVars[i] >= demand[i],
                                    separate = False, modifiable = True))
        newPacking = [0]*len(widths)
        newPacking[i] = numWidthsPerRoll
        packings.append(newPacking)

    # Setting the pricer_data for use in the init and redcost functions
    pricer.data = {}

    pricer.data['cons'] = demandCons

    s.data = {}
    s.data['var'] = PackingVars
    s.data['varnames'] = varNames
    s.data['cons'] = demandCons
    s.data['widths'] = widths
    s.data['demand'] = demand
    s.data['rollLength'] = rollLength
    s.data['packings'] = packings
    s.data['conshdlr'] = conshdlr


    # solve problem
    s.optimize()


    # print original data
    printWidths = '\t'.join(str(e) for e in widths)
    print('\nInput Data')
    print('==========')
    print('Roll Length:', rollLength)
    print('Widths:\t', printWidths)
    print('Demand:\t', '\t'.join(str(e) for e in demand))

    # print solution
    widthOutput = [0]*len(widths)
    print('\nResult')
    print('======')
    print('\t\tSol Value', '\tWidths\t', printWidths)
    for i in range(len(s.data['var'])):
        rollUsage = 0
        solValue = s.getVal(s.data['var'][i]) # round(s.getVal(pricer.data['var'][i]))
        if solValue > 0:
            outline = 'Packing_' + str(i) + ':\t' + str(solValue) + '\t\tCuts:\t '
            for j in range(len(widths)):
                rollUsage += s.data['packings'][i][j]*widths[j]
                widthOutput[j] += s.data['packings'][i][j]*solValue
                outline += str(s.data['packings'][i][j]) + '\t'
            outline += 'Usage:' + str(rollUsage)
            print(outline)

    print('\t\t\tTotal Output:\t', '\t'.join(str(e) for e in widthOutput))

    print(s.getObjVal())


if __name__ == '__main__':
    test_binpacking()

[Steve8881]

I have done more debugging and found something really weird that may caused by the getPyCons() function defined in conshdlr.pxi. When there are more than one branching constraints (i.e., the SameDiff constraints), the constraint list obtained by getPyCons() in the consprop() function is not correct. More interestingly, the result is different from time to time. Sometimes the result is (a traceback object is included)

[<traceback object at 0x7fb363038680>, some_same_name]

Sometimes the result is (a dict is included)

[[{'id': 27, 'name': 'test_binpacking', 'line': 387, 'column': 1, 'source': {'path': '/Users/steve/Desktop/branch-price-and-cut/bin packing.py', 'sourceReference': 0}}], some_same_name]

or (a empty list is included)

[[], some_same_name]

Here is the code:

from pyscipopt import *
import numpy as np
import inspect

entering = lambda: print("IN ", inspect.stack()[1][3]) # debugging magic
leaving = lambda: print("OUT ", inspect.stack()[1][3]) # debugging magic

class SameDiff(Conshdlr):

    def consdataCreate(self, name, item1, item2, constype, node):
        entering()
        cons = self.model.createCons(self, name, stickingatnode = True)
        cons.data = {}
        cons.data['item1'] = item1
        cons.data['item2'] = item2
        cons.data['type'] = constype
        cons.data['npropagatedvars'] = 0
        cons.data['npropagations'] = 0
        cons.data['propagated'] = False
        cons.data['node'] = node
        leaving()
        return cons
    
    def checkVariable(self, cons, var, varid, nfixedvars):
        cutoff = False
        if var.getUbLocal() < 0.5:
            return nfixedvars, cutoff

        constype = cons.data['type']
        packings = self.model.data['packings']
        packingId = varid
        existitem1 = packings[packingId][cons.data['item1']]
        existitem2 = packings[packingId][cons.data['item2']]
        
        if (constype == 'same' and (existitem1 != existitem2)) or (constype == 'diff' and existitem1 and existitem2):
            infeasible, fixed = self.model.fixVar(var, 0.0)
            if infeasible:
                cutoff = True
            else:
                nfixedvars += 1
        return nfixedvars, cutoff

    def consdataFixVariables(self, cons, result):
        nfixedvars = 0
        cutoff = False
        v = cons.data['npropagatedvars']
        while (v < len(self.model.data['var'])) and (not cutoff):
            nfixedvars, cutoff = self.checkVariable(cons, self.model.data['var'][v], v, nfixedvars)
            v += 1
        if (cutoff):
            return {"result": SCIP_RESULT.CUTOFF}
        elif (nfixedvars > 0):
            return {"result": SCIP_RESULT.REDUCEDDOM}
        else:
            return result
        
    def consactive(self, constraint):
        entering()        
        if constraint.data['npropagatedvars'] != len(self.model.data['var']):
            constraint.data['propagated'] = False  
            self.model.repropagateNode(constraint.data['node'])
        leaving()

    def consdeactive(self, constraint):
        entering()
        constraint.data['npropagatedvars'] = len(self.model.data['var'])
        leaving()

    def consprop(self, constraints, nusefulconss, nmarkedconss, proptiming):
        entering()
        result = {"result": SCIP_RESULT.DIDNOTFIND}
        for c in constraints:
            if not c.data['propagated']:
                result = self.consdataFixVariables(c, result)
                c.data['npropagations'] += 1
                if result['result'] != SCIP_RESULT.CUTOFF:
                    c.data['propagated'] = True
                    c.data['npropagatedvars'] = len(self.model.data['var'])
                else:
                    leaving()
                    return {"result": SCIP_RESULT.CUTOFF}
        leaving()
        return result

class MyRyanFosterBranching(Branchrule):

    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):
        entering()
        
        lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars = self.model.getLPBranchCands()

        fractionalities = np.zeros([len(self.model.data['widths']), len(self.model.data['widths'])])

        for ii in range(len(lpcands)):

            solval = lpcandssol[ii]

            packingNum = self.model.data['varnames'].index(lpcands[ii].name)

            currentpacking = self.model.data['packings'][packingNum]

            con_id = np.where(np.array(currentpacking) == 1)[0]

            for i in con_id:
                for j in con_id:
                    if i == j:
                        continue
                    fractionalities[i, j] += solval

        temp = np.fmin(fractionalities, 1 - fractionalities)

        x_index, y_index = np.where(temp == np.max(temp))

        row1 = x_index[0]

        row2 = y_index[0]

        childsame = self.model.createChild(0.0, self.model.getObjVal())

        childdiff = self.model.createChild(0.0, self.model.getObjVal())

        conssame = self.model.data['conshdlr'].consdataCreate("some_same_name", row1, row2, "same", childsame)

        consdiff = self.model.data['conshdlr'].consdataCreate("some_diff_name", row1, row2, "diff", childdiff)

        self.model.addConsNode(childsame, conssame)
        
        self.model.addConsNode(childdiff, consdiff)

        leaving()

        return {"result": SCIP_RESULT.BRANCHED}
        
class MyVarBranching(Branchrule):

    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):

        lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars = self.model.getLPBranchCands()

        integral = lpcands[0]

        self.model.branchVar(integral)

        return {"result": SCIP_RESULT.BRANCHED}
        

class PackingPricer(Pricer):

    # The reduced cost function for the variable pricer
    def pricerredcost(self):

        # Retrieving the dual solutions
        dualSolutions = []
        for i, c in enumerate(self.data['cons']):
            dualSolutions.append(self.model.getDualsolLinear(c))

        #################################################################
        #################################################################

        # Solve the MIP directly

        # Building a MIP to solve the subproblem
        subMIP = Model("BinPacking-Sub")

        # Setting the verbosity level to 0
        subMIP.hideOutput()

        binWidthVars = []
        varNames = []
        varBaseName = "Packing"

        # Variables for the subMIP
        for i in range(len(dualSolutions)):
            varNames.append(varBaseName + "_" + str(i))
            binWidthVars.append(subMIP.addVar(varNames[i], vtype = "B", obj = -1.0 * dualSolutions[i]))

        # Adding the knapsack constraint
        subMIP.addCons(quicksum(w*v for (w,v) in zip(self.model.data['widths'], binWidthVars)) <= self.model.data['rollLength'])

        for ii in range(len(self.model.data['var'])):
            if self.model.data['var'][ii].getUbLocal() < 0.5:
                currentPacking = self.model.data['packings'][ii]
                subMIP.addCons(quicksum(w*v for (w,v) in zip(currentPacking, binWidthVars)) <= np.sum(currentPacking) - 1)


        # Solving the subMIP to generate the most negative reduced cost packing
        subMIP.optimize()

        objval = 1 + subMIP.getObjVal()

        #################################################################
        #################################################################
        

        # Adding the column to the master problem
        if objval < -1e-08:
            currentNumVar = len(self.model.data['var'])

            # Creating new var; must set pricedVar to True
            newVar = self.model.addVar("NewPacking_" + str(currentNumVar), vtype = "I", obj = 1.0, pricedVar = True)

            self.model.data['varnames'].append("NewPacking_" + str(currentNumVar))

            # Adding the new variable to the constraints of the master problem
            newPacking = []
            for i, c in enumerate(self.data['cons']):
                coeff = round(subMIP.getVal(binWidthVars[i]))
                self.model.addConsCoeff(c, newVar, coeff)

                newPacking.append(coeff)
            
            self.model.data['packings'].append(newPacking)
            self.model.data['var'].append(newVar)

        return {'result':SCIP_RESULT.SUCCESS}

    # The initialisation function for the variable pricer to retrieve the transformed constraints of the problem
    def pricerinit(self):
        for i, c in enumerate(self.data['cons']):
            self.data['cons'][i] = self.model.getTransformedCons(c)


def test_binpacking():

    s = Model()

    s.setPresolve(SCIP_PARAMSETTING.OFF)

    s.setSeparating(SCIP_PARAMSETTING.OFF)

    s.setHeuristics(SCIP_PARAMSETTING.OFF)


    pricer = PackingPricer(priority = 5000000)
    s.includePricer(pricer, "BinPackingPricer", "Pricer")

    conshdlr = SameDiff()

    s.includeConshdlr(conshdlr, "SameDiff", "SameDiff Constraint Handler",
                          chckpriority=2000000, propfreq=1)

    my_branchrule = MyRyanFosterBranching(s)

    # my_branchrule = MyVarBranching(s)

    s.includeBranchrule(my_branchrule, "test branch", "test branching and probing and lp functions",
                    priority=10000000, maxdepth= -1, maxbounddist=1)

    # item widths
    widths = [71,80,76,83,41,78,70,23,42,87,43,84,60,55,49]
    # width demand
    demand = [1] * 15 

    # roll length
    rollLength = 150
    assert len(widths) == len(demand)

    # adding the initial variables
    cutpackingVars = []
    varNames = []
    varBaseName = "packing"
    packings = []

    for i in range(len(widths)):
        varNames.append(varBaseName + "_" + str(i))
        cutpackingVars.append(s.addVar(varNames[i], vtype = 'B', obj = 1.0))

    demandCons = []
    for i in range(len(widths)):

        numWidthsPerRoll = 1 
        demandCons.append(s.addCons(numWidthsPerRoll*cutpackingVars[i] >= demand[i],
                                    separate = False, modifiable = True))
        newpacking = [0]*len(widths)
        newpacking[i] = numWidthsPerRoll
        packings.append(newpacking)

    pricer.data = {}

    pricer.data['cons'] = demandCons

    s.data = {}
    s.data['var'] = cutpackingVars
    s.data['varnames'] = varNames
    s.data['cons'] = demandCons
    s.data['widths'] = widths
    s.data['demand'] = demand
    s.data['rollLength'] = rollLength
    s.data['packings'] = packings
    s.data['conshdlr'] = conshdlr

    s.optimize()

    # print original data
    printWidths = '\t'.join(str(e) for e in widths)
    print('\nInput Data')
    print('==========')
    print('Roll Length:', rollLength)
    print('Widths:\t', printWidths)
    print('Demand:\t', '\t'.join(str(e) for e in demand))

    # print solution
    widthOutput = [0]*len(widths)
    print('\nResult')
    print('======')
    print('\t\tSol Value', '\tWidths\t', printWidths)
    for i in range(len(s.data['var'])):
        rollUsage = 0
        solValue = s.getVal(s.data['var'][i]) 
        if solValue > 0:
            outline = 'Packing_' + str(i) + ':\t' + str(solValue) + '\t\tCuts:\t '
            for j in range(len(widths)):
                rollUsage += s.data['packings'][i][j]*widths[j]
                widthOutput[j] += s.data['packings'][i][j]*solValue
                outline += str(s.data['packings'][i][j]) + '\t'
            outline += 'Usage:' + str(rollUsage)
            print(outline)

    print('\t\t\tTotal Output:\t', '\t'.join(str(e) for e in widthOutput))

    print(s.getObjVal())


if __name__ == '__main__':
    test_binpacking()

[fserra]

Ok I think I know what is going on. I checked the code on https://github.com/SCIP-Interfaces/PySCIPOpt/blob/use-probdata-to-store-model/tests/test_coloring.py and there we had a function addPyConsNode which is basically the same as addConsNode, except for a Py_INCREF(cons) missing in addConsNode.

After adding that line it fixes the problem here. @purin0604 can you check whether that also fixes your issue?

The problem is that, with a python constraint handler, one calls Model.createCons. This creates a constraint in SCIP and a python Constraint object. Then, one uses the python Constraint object to add the constraint to the node.
This will usually happen in the scope of a callback (like branchexeclp in the MyRyanFosterBranching code above), so after this function is called the python Constraint object is not referenced by no one inside python. It is however, reference by SCIP, because SCIP stores a pointer to the python Constraint object in consdata of the SCIP constraint created in Model.createCons.
So then this C pointer to the python object is transformed back to a python Constraint getPyCons, but it points probably to garbage, which generates the segfault.

So by increasing the reference to the python Constraint object, it should prevent python from deleting it. The reference is later decreased in the delete callback (see PyConsDelete in conshdlr.pxi)

@mattmilten should we rename the function addConsNode and addConsLocal to addPyConsNode and addPyConsLocal?
This would be analogous to addPyCons.

@Steve8881 are you trying to translate the BPP example from SCIP? In any case, would you like to share the example with us? There has been some demand to have a more complete branch and price example, but I haven't find the time to adapt the coloring one nor to translate the BPP from SCIP.

[fserra]

btw, after this other problems occur in your code @Steve8881. The callback consprop in your implementation is missing the return of a dictionary specifying the result

[mattmilten]

I am fine with renaming addConsNode and addConsLocal. This might introduce issues with older code but PySCIPOpt never claimed to be backward compatible, anyway.

[Steve8881]

Hi @fserra, thanks for your reply! Of course, I would like to share my code with you but there are still some small problems need to solve.

After adding Py_INCREF(cons) to the addConsNode() function, and fixing some typos, the algorithm works. However, after the branch-and-price process and obtaining the solution, an assertion is failed: Assertion failed: ((*conshdlr)->nconss == 0), function SCIPconshdlrFree, file /Users/steve/scipoptsuite-7.0.1/scip/src/scip/cons.c, line 2351. Abort trap: 6.

Besides, I find that in the SCIP example, the conscheck() method of the SameDiff constraint handler is set as NULL but this method is called frequently in my PySCIPOpt implementation. When I try to do something, for example, checking whether a solution is feasible or not in the conscheck() method, another assertion is failed: Assertion failed: (scip != NULL), function SCIPcheckStage, file /Users/steve/scipoptsuite-7.0.1/scip/src/scip/debug.c, line 2044. Abort trap: 6. However, when the conscheck() method just returns {"result": SCIP_RESULT.FEASIBLE} directly, everything works well. Even the result of the algorithm is correct. What is going on here?

Finally, there are also some questions about the SCIP and the branch-and-price algorithm for BPP.

• In the pricing problem of the BPP example, besides the knapsack constraint, the branching decision constraints and the fixed variable constraints (avoid to generate columns which are fixed to zero) are added to the sub-MIP. Note that the variables are fixed to zero according to the branching decisions, don't the latter two types of constraints mean the same thing? However, when I remove the fixed variable constraints from the pricing problem, the running time of the algorithm increases sharply but the solution is still correct, which confuses me a lot. Why?

• What is the default node selector in SCIP? In the BPP example, the priorities of all nodes are set as 0.0. What happens in this case?

• If I want to use one of the node selector implemented by SCIP, e.g., dfs, is there some function to give the child node a priority? What if a custom node selector is applied (especially when the model is not branched by a variable but constraints)?

Thank you.

Steve

PS: Here is the code

from pyscipopt import *
import numpy as np
import inspect


def entering():
    return print("IN ", inspect.stack()[1][3])  # debugging magic


def leaving():
    return print("OUT ", inspect.stack()[1][3])  # debugging magic


class SameDiff(Conshdlr):
    def consdataCreate(self, name, item1, item2, constype, node):
        cons = self.model.createCons(self, name, stickingatnode=True)
        cons.data = {}
        cons.data["item1"] = item1
        cons.data["item2"] = item2
        cons.data["type"] = constype
        cons.data["npropagatedvars"] = 0
        cons.data["npropagations"] = 0
        cons.data["propagated"] = False
        cons.data["node"] = node
        return cons

    def checkVariable(self, cons, var, varid, nfixedvars):
        cutoff = False
        if var.getUbLocal() < 0.5:
            return nfixedvars, cutoff

        constype = cons.data["type"]
        packings = self.model.data["packings"]
        packingId = varid
        existitem1 = packings[packingId][cons.data["item1"]]
        existitem2 = packings[packingId][cons.data["item2"]]

        if (constype == "same" and (existitem1 != existitem2)) or (constype == "diff" and existitem1 and existitem2):
            infeasible, fixed = self.model.fixVar(var, 0.0)
            if infeasible:
                cutoff = True
            else:
                nfixedvars += 1
        return nfixedvars, cutoff

    def consdataFixVariables(self, cons, result):
        nfixedvars = 0
        cutoff = False
        v = cons.data["npropagatedvars"]
        while (v < len(self.model.data["var"])) and (not cutoff):
            nfixedvars, cutoff = self.checkVariable(cons, self.model.data["var"][v], v, nfixedvars)
            v += 1
        if cutoff:
            return {"result": SCIP_RESULT.CUTOFF}
        elif nfixedvars > 0:
            return {"result": SCIP_RESULT.REDUCEDDOM}
        else:
            return result

    def consactive(self, constraint):
        # entering()
        if constraint.data["npropagatedvars"] != len(self.model.data["var"]):
            constraint.data["propagated"] = False
            self.model.repropagateNode(constraint.data["node"])
        # leaving()

    def consdeactive(self, constraint):
        # entering()
        constraint.data["npropagatedvars"] = len(self.model.data["var"])
        # leaving()

    def conscheck(
        self,
        constraints,
        solution,
        checkintegrality,
        checklprows,
        printreason,
        completely,
    ):
        # for cons in constraints:
        #     item1 = cons.data["item1"]
        #     item2 = cons.data["item1"]
        #     packings = self.model.data["packings"]
        #     for i in range(len(self.model.data["var"])):
        #         sol = solution[self.model.data["var"][i]]
        #         if sol < 0.5:
        #             continue
        #         else:
        #             packingId = i
        #             existitem1 = packings[packingId][item1]
        #             existitem2 = packings[packingId][item2]
        #             if cons.data['type'] == 'same' and existitem1 != existitem2:
        #                 return {"result": SCIP_RESULT.INFEASIBLE}
        #             elif cons.data['type'] == 'diff' and existitem1 and existitem2:
        #                 return {"result": SCIP_RESULT.INFEASIBLE}

        return {"result": SCIP_RESULT.FEASIBLE}

    def consenfolp(self, constraints, nusefulconss, solinfeasible):
        pass

    def consenfops(self, constraints, nusefulconss, solinfeasible, objinfeasible):
        pass

    def consprop(self, constraints, nusefulconss, nmarkedconss, proptiming):
        result = {"result": SCIP_RESULT.DIDNOTFIND}
        for c in constraints:
            if not c.data["propagated"]:
                result = self.consdataFixVariables(c, result)
                c.data["npropagations"] += 1
                if result["result"] != SCIP_RESULT.CUTOFF:
                    c.data["propagated"] = True
                    c.data["npropagatedvars"] = len(self.model.data["var"])
                else:
                    leaving()
                    return {"result": SCIP_RESULT.CUTOFF}
        return result


class MyRyanFosterBranching(Branchrule):
    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):
        lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars = self.model.getLPBranchCands()

        fractionalities = np.zeros([len(self.model.data["widths"]), len(self.model.data["widths"])])

        for ii in range(len(lpcands)):

            solval = lpcandssol[ii]

            packingNum = self.model.data["varnames"].index(lpcands[ii].name)

            currentpacking = self.model.data["packings"][packingNum]

            con_id = np.where(np.array(currentpacking) == 1)[0]

            for i in con_id:
                for j in con_id:
                    if i == j:
                        continue
                    fractionalities[i, j] += solval

        temp = np.fmin(fractionalities, 1 - fractionalities)

        x_index, y_index = np.where(temp == np.max(temp))

        row1 = x_index[0]

        row2 = y_index[0]

        childsame = self.model.createChild(0.0, self.model.getLocalTransEstimate())

        childdiff = self.model.createChild(0.0, self.model.getLocalTransEstimate())

        conssame = self.model.data["conshdlr"].consdataCreate("some_same_name", row1, row2, "same", childsame)

        consdiff = self.model.data["conshdlr"].consdataCreate("some_diff_name", row1, row2, "diff", childdiff)

        self.model.addConsNode(childsame, conssame)

        self.model.addConsNode(childdiff, consdiff)

        self.model.data['branchingCons'].append(conssame)

        self.model.data['branchingCons'].append(consdiff)

        return {"result": SCIP_RESULT.BRANCHED}


class MyVarBranching(Branchrule):
    def __init__(self, model):
        self.model = model

    def branchexeclp(self, allowaddcons):

        (
            lpcands,
            lpcandssol,
            lpcandsfrac,
            nlpcands,
            npriolpcands,
            nfracimplvars,
        ) = self.model.getLPBranchCands()

        integral = lpcands[0]

        self.model.branchVar(integral)

        return {"result": SCIP_RESULT.BRANCHED}


class Pricer(Pricer):
    def addBranchingDecisionConss(self, subMIP, binWidthVars):
        for cons in self.model.data['branchingCons']:
            if (not cons.isActive()):
                continue
            item1 = cons.data['item1']
            item2 = cons.data['item2']

            if cons.data['type'] == 'same':
                subMIP.addCons(binWidthVars[item1] - binWidthVars[item2] == 0)
            elif cons.data['type'] == 'diff':
                subMIP.addCons(binWidthVars[item1] + binWidthVars[item2] <= 1)

    def addFixedVarsConss(self, subMIP, binWidthVars):
        for ii in range(len(self.model.data["var"])):
            if self.model.data["var"][ii].getUbLocal() < 0.5:
                currentPacking = self.model.data["packings"][ii]
                subMIP.addCons(
                    quicksum(w * v for (w, v) in zip(currentPacking, binWidthVars)) <= np.sum(currentPacking) - 1)

    # The reduced cost function for the variable pricer
    def pricerredcost(self):

        # Retrieving the dual solutions
        dualSolutions = []
        for i, c in enumerate(self.data["cons"]):
            dualSolutions.append(self.model.getDualsolLinear(c))
            # dualSolutions.append(self.model.getDualsolSetppc(c))

        #################################################################
        #################################################################

        # Solve the MIP directly

        # Building a MIP to solve the subproblem
        subMIP = Model("BinPacking-Sub")

        # Setting the verbosity level to 0
        subMIP.hideOutput()

        binWidthVars = []
        varNames = []
        varBaseName = "Packing"

        # Variables for the subMIP
        for i in range(len(dualSolutions)):
            varNames.append(varBaseName + "_" + str(i))
            binWidthVars.append(subMIP.addVar(varNames[i], vtype="B", obj=-1.0 * dualSolutions[i]))

        # Adding the knapsack constraint
        subMIP.addCons(
            quicksum(w * v for (w, v) in zip(self.model.data["widths"], binWidthVars)) <= self.model.data["rollLength"])

        # Avoid to generate columns which are fixed to zero
        self.addFixedVarsConss(subMIP, binWidthVars)

        # Add branching decisions constraints to the sub SCIP
        self.addBranchingDecisionConss(subMIP, binWidthVars)

        # Solving the subMIP to generate the most negative reduced cost packing
        subMIP.optimize()

        objval = 1 + subMIP.getObjVal()

        #################################################################
        #################################################################

        # Adding the column to the master problem
        if objval < -1e-08:

            currentNumVar = len(self.model.data["var"])

            # Creating new var; must set pricedVar to True
            newVar = self.model.addVar("NewPacking_" + str(currentNumVar), vtype="I", obj=1.0, pricedVar=True)

            self.model.data["varnames"].append("NewPacking_" + str(currentNumVar))

            # Adding the new variable to the constraints of the master problem
            newPacking = []
            for i, c in enumerate(self.data["cons"]):
                coeff = round(subMIP.getVal(binWidthVars[i]))
                self.model.addConsCoeff(c, newVar, coeff)
                # self.model.addVarBasicSetcover(c, newVar)

                newPacking.append(coeff)

            self.model.data["packings"].append(newPacking)
            self.model.data["var"].append(newVar)

        return {"result": SCIP_RESULT.SUCCESS}

    # The initialisation function for the variable pricer to retrieve the transformed constraints of the problem
    def pricerinit(self):
        for i, c in enumerate(self.data["cons"]):
            self.data["cons"][i] = self.model.getTransformedCons(c)


def test_binpacking():

    s = Model()

    # s.setPresolve(SCIP_PARAMSETTING.OFF)

    s.setIntParam("presolving/maxrestarts", 0)

    s.setParam('limits/time', 30 * 60)

    s.setSeparating(SCIP_PARAMSETTING.OFF)

    # s.setHeuristics(SCIP_PARAMSETTING.OFF)

    s.setObjIntegral()

    pricer = Pricer(priority=5000000)
    s.includePricer(pricer, "BinPackingPricer", "Pricer")

    conshdlr = SameDiff()

    s.includeConshdlr(
        conshdlr,
        "SameDiff",
        "SameDiff Constraint Handler",
        chckpriority=2000000,
        propfreq=1,
    )

    my_branchrule = MyRyanFosterBranching(s)

    # my_branchrule = MyVarBranching(s)

    s.includeBranchrule(
        my_branchrule,
        "test branch",
        "test branching and probing and lp functions",
        priority=10000000,
        maxdepth=-1,
        maxbounddist=1,
    )

    # item widths
    widths = [
        42, 69, 67, 57, 93, 90, 38, 36, 45, 42, 33, 79, 27, 57, 44, 84, 86, 92, 46, 38, 85, 33, 82, 73, 49, 70, 59, 23,
        57, 72, 74, 69, 33, 42, 28, 46, 30, 64, 29, 74, 41, 49, 55, 98, 80, 32, 25, 38, 82, 30, 35, 39, 57, 84, 62, 50,
        55, 27, 30, 36, 20, 78, 47, 26, 45, 41, 58, 98, 91, 96, 73, 84, 37, 93, 91, 43, 73, 85, 81, 79, 71, 80, 76, 83,
        41, 78, 70, 23, 42, 87, 43, 84, 60, 55, 49, 78, 73, 62, 36, 44, 94, 69, 32, 96, 70, 84, 58, 78, 25, 80, 58, 66,
        83, 24, 98, 60, 42, 43, 43, 39
    ]
    # width demand
    demand = [1] * 120

    # roll length
    rollLength = 150
    assert len(widths) == len(demand)

    # adding the initial variables
    binpackingVars = []
    varNames = []
    varBaseName = "Packing"
    packings = []

    for i in range(len(widths)):
        varNames.append(varBaseName + "_" + str(i))
        binpackingVars.append(s.addVar(varNames[i], vtype="B", obj=1.0))

    demandCons = []
    for i in range(len(widths)):

        numWidthsPerRoll = 1
        demandCons.append(
            s.addCons(
                numWidthsPerRoll * binpackingVars[i] >= demand[i],
                separate=False,
                modifiable=True,
            ))
        newpacking = [0] * len(widths)
        newpacking[i] = numWidthsPerRoll
        packings.append(newpacking)

    pricer.data = {}

    pricer.data["cons"] = demandCons

    s.data = {}
    s.data["var"] = binpackingVars
    s.data["varnames"] = varNames
    s.data["cons"] = demandCons
    s.data["widths"] = widths
    s.data["demand"] = demand
    s.data["rollLength"] = rollLength
    s.data["packings"] = packings
    s.data["conshdlr"] = conshdlr
    s.data['branchingCons'] = []

    s.optimize()

    # print original data
    printWidths = "\t".join(str(e) for e in widths)
    print("\nInput Data")
    print("==========")
    print("Roll Length:", rollLength)
    print("Widths:\t", printWidths)
    print("Demand:\t", "\t".join(str(e) for e in demand))

    # print solution
    widthOutput = [0] * len(widths)
    print("\nResult")
    print("======")
    print("\t\tSol Value", "\tWidths\t", printWidths)
    for i in range(len(s.data["var"])):
        rollUsage = 0
        solValue = s.getVal(s.data["var"][i])
        if solValue > 0:
            outline = "Packing_" + str(i) + ":\t" + \
                str(solValue) + "\t\tCuts:\t "
            for j in range(len(widths)):
                rollUsage += s.data["packings"][i][j] * widths[j]
                widthOutput[j] += s.data["packings"][i][j] * solValue
                outline += str(s.data["packings"][i][j]) + "\t"
            outline += "Usage:" + str(rollUsage)
            print(outline)

    print("\t\t\tTotal Output:\t", "\t".join(str(e) for e in widthOutput))

    print(s.getObjVal())


if __name__ == "__main__":
    test_binpacking()

[fserra]

Hi @Steve8881

the problem with the assert at the end is an old one (see #26). It is a technicality, but the issue here is that the constraints generated by your constraint handler are created and then added, but nobody releases it (in contrast the constraints generated by SCIP's constraint handler are created and added in the same function call and we release them once in there). I am still not 100% sure how to solve this issue.

The problem with conscheck is that we, in python we have dummy implementations of all callbacks. We can't know a priori which callback is useful or not, so we need them all. Thus, conscheck is going to get called because it is not null.
The scip != NULL assert is a nice one. You are calling solution[self.model.data["var"][i]], which uses the SCIP stored inside the solution object. This solution object is created in conshdlr.pxi

cdef SCIP_RETCODE PyConsCheck (SCIP* scip, SCIP_CONSHDLR* conshdlr, SCIP_CONS** conss, int nconss, SCIP_SOL* sol, SCIP_Bool checkintegrality,
                               SCIP_Bool checklprows, SCIP_Bool printreason, SCIP_Bool completely, SCIP_RESULT* result):
    PyConshdlr = getPyConshdlr(conshdlr)
    cdef constraints = []
    for i in range(nconss):
        constraints.append(getPyCons(conss[i]))
    solution = Solution()
    solution.sol = sol
    result_dict = PyConshdlr.conscheck(constraints, solution, checkintegrality, checklprows, printreason, completely)
    result[0] = result_dict.get("result", <SCIP_RESULT>result[0])
    return SCIP_OKAY

However, we are not setting the scip inside the solution, hence the problem. The fix is just to add solution.scip = scip. I will fix this.

The rest of your questions I cannot answer properly. Maybe @stephenjmaher knows more. Just some guesses:

I would imagine that there are more fixings apart from the branching decisions.

The default node selector is estimate. When comparing two nodes, it checks their estimate. if the estimate are equal then it checks the lower bound of the nodes, if this is equal then it does some check depending on the type of the node and depth. However, if everything is equal, the comparison returns 0. I am not sure how are they sorten then, if randomly or not.

when you create a child you can give it a priority.

[purin0604]

@fserra Thank you!! Adding Py_INCREF(cons) solved my problem as well.

I also have problem python3: /opt/scipoptsuite-7.0.0/scip/src/scip/cons.c:2351: SCIPconshdlrFree: Assertion (*conshdlr)->nconss == 0' failed.` as Steve8881.
If this is just a technical issue, can I ignore the assertion error and be assured with the optimal solution?

[stephenjmaher]

Just to follow up on the questions from @Steve8881, my understanding is the same as @fserra.

In regards to the additional constraints in the pricer. The constraints setting the variables to 0 are to capture any other variable fixings that have occurred that are not from branching decisions. There are a lot of different ways that the upper bound of a variable can be reduced, so this is trying to incorporate some of that into the pricer.

[fserra]

@purin0604, yes, that doesn't affect the correctness of the algorithm

[Steve8881]

Hi @fserra and @stephenjmaher

Thanks for your very quick replies! Adding solution.scip = scip fixes my problem.

I do know that I could give the child node a priority when creating it but the value of the priority is not easy to determine. The value is related to the specific node selector we use. For example, if depth-first-search (DFS) is used, we should give deeper nodes higher priorities when creating the child nodes in branchexeclp() function. It follows that we should know the position of the current node in the searching tree, which is not easy to obtain. So I guess there may be some method defined in the node selector class to calculate the priorities we need to give the nodes.

In your description of the estimate node selector, the priorities of the nodes seem to be useless. Is this true?

Thank you.

Steve

[fserra]

Hi. What do you mean by the position of the current node in the search tree? There is a function to get the depth of a node (getDepth).

As I said before I don't really know how the node selection in SCIP works. But from some reading here https://www.scipopt.org/doc/html/NODESEL.php I understand that the priority given when creating a node influences for example, the function SCIPgetPrioChild, which is used inside the select callback of a node selector. Also, the node selector plugin has its own method for comparing nodes (the compare callback). For example, best child according to the node selector's own criteria can be obtained via SCIPgetBestChild

What I described above was the comparing mechanism of the estimate node selector. Indeed, there the priorities given when the nodes were created seem to be useless. However, when actually selecting a node, this priority is taken into account (i.e. the estimate node selector calls SCIPgetPrioChild)

The strategy for selecting a node of the estimate nodeselector seems to be a mixture between plunging, breadth first search and best bound (you can have a look in nodesel_estimate.c, the code is not so bad)