Adding Lazy Constraints with condition in Gurobi via Python 2.7

[Taner Cokyasar]

Hi everyone,

I am trying to solve a problem by adding cut constraints. I generated a logic to add the constraints, however, I do not know what the most efficient way of adding these cuts is. I read something about lazy constraints and it looks like my constraints can be considered as lazy constraints. I tried to explicitly show my problem in the following link:

https://stackoverflow.com/questions/47383095/adding-lazy-constraints-with-condition-in-gurobi-via-python-2-7

The link includes my code and the logic of constraints addition. If there is a function in Gurobi that can help me to add these constraints in a sequential way, I would be very happy learning it.

Thanks for your help in advance.

Best,

Taner

[Daniel Espinoza]

Hi Taner,

I would suggest to try your constraints as simple user cuts (via a user cut callback), not as lazy, and see if the performance is better. Or, if you have a fixed list of constraints, maybe add them at the beginning of the model and treat them as lazy type 3.

Best,

Daniel

[Taner Cokyasar]

Hi Daniel,

Is there any example code for Python that shows how to use simple user cuts via cut callback?  I have a fixed list of constraints. However, if the condition is not met, I expand the numbers (increase the increments from 0.1 to 0.01 or from 0.01 to 0.0001 in the range 0,1)on my list and solve it again until it meets the condition. In this case, any example that I can follow for lazy type 3?

Best,

Taner

[Daniel Espinoza]

Hi Taner,

You can define the constraint as regular ones, and set the attribute Lazy=3 for them.

The documentation is here

http://www.gurobi.com/documentation/7.5/refman/lazy.html#attr:Lazy

It's a constraint attribute, so you would only set it for those constraints.

For cut callbacks, you can take a look here:

http://www.gurobi.com/documentation/7.5/examples/callback_py.html

For cut callbacks you can query the current relaxation solution and check your list of constraints and add the most violated constraints.

In this case you should not add the constraints from the beginning.

It's not clear which of the two approaches would be better, you should test it.

Best,

Daniel

[Taner Cokyasar]

Dr. Espinoza,

Firstly, thanks for your response. I read both of them and attempted to define a callback function that should supposedly do what I am looking for. However, it is not working. I added my initial constraints in the base model and expected the callback to increase the number of variables and constraints with the list rhoset. The error I get is: Exception gurobipy.GurobiError: GurobiError() in 'gurobipy.callbackstub' ignored

If you may help me to see what I am doing wrong, I would really appreciate it.

Best,

Taner

My code is as follows:

def mycallback(m, where):

    if where == GRB.Callback.MIP:

        # General MIP callback

        objbst = m.cbGet(GRB.Callback.MIP_OBJBST)

        

        lambdaforlisofA = quicksum(lambdaa[j-1]*yA[j].x for j in range(1, J+1)).getValue()

        rhoforlistofA = quicksum((lambdaa[j-1]/mu[j-1])*yA[j].x for j in range(1, J+1)).getValue()

        thetaforlistofA = quicksum((lambdaa[j-1]/(mu[j-1]**2))*yA[j].x for j in range(1, J+1)).getValue()

        WqforlistofA = thetaforlistofA/(2*(1-rhoforlistofA))

        CostofA = quicksum(lambdaa[j-1]*(CA[j-1] + d[j-1]*(1/mu[j-1]+WqforlistofA))*yA[j].x for j in range(1, J+1)).getValue()

        CostofI = quicksum((CI[j-1]*lambdaa[j-1]+h[j-1]*(Sstar[j-1]-lambdaa[j-1]*tauI[j-1])

               + (h[j-1]+d[j-1])*(lambdaa[j-1]*tauI[j-1]*PSstarminusone[j-1]

                  - Sstar[j-1]*PSstar[j-1]))*yI[j].x for j in range(1, J+1)).getValue()

        CostofO = quicksum(((CO[j-1] +d[j-1]*tauO[j-1])*lambdaa[j-1])*yO[j].x for j in range(1, J+1)).getValue()

        TotalCost = CostofA + CostofI + CostofO + f

        if 1-(objbst / TotalCost) > 0.001:

            inc2 = 0.01

            rhoset = sorted(rhoset.append(([round(i,3) for i in list(np.arange(m.getVarByName("rho").x-0.1, m.getVarByName("rho").x+0.1, inc2))])))

            print "The solution gap is currently %s percentage." % (1-(objbst / TotalCost))

        else:    

            print "ObjVal:", m.objVal

            print "rho:", m.getVarByName("rho").x

            print "theta:", m.getVarByName("theta").x

            print "lambda:", m.getVarByName("lambdatotal").x

            print "Wq:", m.getVarByName("Wq").x

            print "x:", m.getVarByName("x").x

            print "The solution gap is %s percentage." % (1-(objbst / TotalCost))

            if 1-(objbst / TotalCost) > 0.001:

                inc2 = 0.01

                rhoset = sorted(rhoset.append(([round(i,3) for i in list((np.arange(m.getVarByName("rho").x-0.1, m.getVarByName("rho").x+0.1, inc2)))])))

                print "The solution gap is currently %s percentage." % (1-(objbst / TotalCost))

            else:

                print "ObjVal:", m.objVal

                print "rho:", m.getVarByName("rho").x

                print "theta:", m.getVarByName("theta").x

                print "lambda:", m.getVarByName("lambdatotal").x

                print "Wq:", m.getVarByName("Wq").x

                print "x:", m.getVarByName("x").x

                print "The solution gap is %s percentage." % (1-(objbst / TotalCost))

                if 1-(objbst / TotalCost) > 0.001:

                    inc3 = 0.001

                    rhoset = sorted(rhoset.append(([round(i,3) for i in list((np.arange(m.getVarByName("rho").x, m.getVarByName("rho").x, inc3)))])))

                    print "The solution gap is currently %s percentage." % (1-(objbst / TotalCost))

                else:    

                    print "ObjVal:", m.objVal

                    print "rho:", m.getVarByName("rho").x

                    print "theta:", m.getVarByName("theta").x

                    print "lambda:", m.getVarByName("lambdatotal").x

                    print "Wq:", m.getVarByName("Wq").x

                    print "x:", m.getVarByName("x").x

                    print "The solution gap is %s percentage." % (1-(objbst / TotalCost))

############################################# Create the Model ############################################

m = Model("ExampleCode")

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

############################################## Add Variables ##############################################

inc = 0.1

lowerinc = 0

upperinc = 1

rhoset = sorted([round(i,15) for i in list((np.arange(lowerinc, upperinc, inc)))] + [0.99])

delta = {}  #THESE ARE THE SET OF VARIABLES THAT WILL BE ADDED IF THE SOLUTION IS NOT AS DESIRED

for i in range(1, len(rhoset)+1):

    delta[i] = m.addVar(vtype=GRB.BINARY, name = "delta%s" % str([i]))

yA = {}

yI = {}

yO = {}

alpha = {}

x = m.addVar(vtype=GRB.BINARY, name = "x")

lambdatotal = m.addVar(vtype=GRB.CONTINUOUS, name = "lambdatotal")

rho = m.addVar(vtype=GRB.CONTINUOUS, name = "rho")

theta = m.addVar(vtype=GRB.CONTINUOUS, name = "theta")

Wq = m.addVar(vtype=GRB.CONTINUOUS, name = "Wq")

for j in range(1, J+1):

    yA[j] = m.addVar(vtype=GRB.BINARY, name = "yA%s" % str([j]))

    yI[j] = m.addVar(vtype=GRB.BINARY, name = "yI%s" % str([j]))

    yO[j] = m.addVar(vtype=GRB.BINARY, name = "yO%s" % str([j]))

    alpha[j] = m.addVar(vtype=GRB.CONTINUOUS, name = "alpha%s" % str([j]))

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

############################################## Add Objective ##############################################

m.setObjective(f*x + quicksum(lambdaa[j-1]*(CA[j-1] + d[j-1]*(1/mu[j-1]))*yA[j] for j in range(1, J+1))

   + quicksum(alpha[j] for j in range(1, J+1))

     + quicksum((CI[j-1]*lambdaa[j-1]+h[j-1]*(Sstar[j-1]-lambdaa[j-1]*tauI[j-1])

       + (h[j-1]+d[j-1])*(lambdaa[j-1]*tauI[j-1]*PSstarminusone[j-1]

         - Sstar[j-1]*PSstar[j-1]))*yI[j] for j in range(1, J+1))

           + quicksum(((CO[j-1]+d[j-1]*tauO[j-1])*lambdaa[j-1])*yO[j] for j in range(1, J+1)), GRB.MINIMIZE)

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

############################################# Add Constraints #############################################

for j in range(1, J+1):

    m.addConstr(yA[j] + yI[j] + yO[j] == 1, name = "Cons2")                                                                                 #Constraint Set 2

m.addConstr(quicksum(yA[j] for j in range(1, J+1)) <= J*x, name = "Cons3")                                                     #Constraint 3

m.addConstr(quicksum(lambdaa[j-1]*yA[j] for j in range(1, J+1)) == lambdatotal, name = "Cons4")                  #Constraint 4

m.addConstr(quicksum(lambdaa[j-1]/mu[j-1]*yA[j] for j in range(1, J+1)) == rho, name = "Cons5")                   #Constraint 5

m.addConstr(quicksum(lambdaa[j-1]/mu[j-1]/mu[j-1]*yA[j] for j in range(1, J+1)) == theta, name = "Cons6")    #Constraint 6

for j in range(1, J+1):

    m.addConstr(alpha[j] >= Wq*d[j-1]*lambdaa[j-1]-(1-yA[j])*M1, name = "Cons7")                                           #Constraint 7

m.addConstr(rho <= 0.99, name = "Cons8")                                                                                                       #Constraint 8

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

########################################### Additional Constraints ###########################################

for i in range(len(rhoset)): #THESE ARE THE SET OF CONSTRAINTS WILL BE ADDED IF THE SOLUTION IS NOT AS DESIRED

    m.addConstr(rho - rhoset[i] + epsilon <= delta[i+1], name = "Cons9")                                                           #Constraint 9

    m.addConstr(Wq >= (1/(2*(1-rhoset[i])))*theta - M2*(1-delta[i+1]), name = "Cons10")                                  #Constraint 10

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

m.update()

m.optimize(mycallback)

print "ObjVal:", m.objVal

print "rho:", m.getVarByName("rho").x

print "theta:", m.getVarByName("theta").x

print "lambda:", m.getVarByName("lambdatotal").x

print "Wq:", m.getVarByName("Wq").x

print "x:", m.getVarByName("x").x

[Daniel Espinoza]

Hi Taner,

You can call me Daniel.

I did not read your code in detail, however, there seems to be a misunderstanding.

First: During a callback you can not add variables to the model.

Second: it seems to me that you could add all constraints from the beginning, without callback or lazy attribute needed and solve a monolithic model. Specially so because the number of extra constraints don't seem to be large.

At least it would be worth it a try. If what you are trying to do is to implement a custom stop criteria, that could be done within a callback.

Finally, as a recommendation, it's much more efficient to reference variables by objects than to query them by name. Remember that addVars and addConstrs will return a tupledict of variables and constraints respectively.

Best,

Daniel