Skip to main content

Differentiating an equation and updating as a constraint

Answered
Varun Deshpande
  • Gurobi-versary
  • First Question
  • First Comment

How should an expression be partially differentiated and added as a constraint to the model? Clearly using sympy to differntiate isn't too helpful since.addConstr() has not accepted the constraint. The solution of the problem below should be 0 and not infeasible/unbounded.

m = gp.Model()
x = m.addVar(name = 'x', lb = 0, ub = 3)
y = m.addVar(name = 'y', lb = -gp.GRB.INFINITY, ub = gp.GRB.INFINITY)

# Enter constraints and objective functions
constr_lower = [x-y>=1, y <= 0, x + y <= 2]
constr_upper = []
F_upper = (x-1.5)**2 + (y+1)**2
f_lower = -y

for cs in constr_lower:
  m.addConstr(cs)

mu_1 = sym.Symbol('mu_1')
mu_2 = sym.Symbol('mu_2')
lambda_1 = sym.Symbol('lambda_1')
lambda_2 = sym.Symbol('lambda_2')
beta = sym.Symbol('beta')
gamma = sym.Symbol('gamma')
lagrangian = mu_1*phi[z[0]] +mu_2*phi[z[1]]-lambda_1*mu_1-lambda_2*mu_2 + beta *(mu_1+mu_2 - 1) + gamma*(mu_1*z[0] + mu_2*z[1])

c1 = sym.diff(lagrangian, mu_1)
c2 = sym.diff(lagrangian, mu_2)
c3 = sym.diff(lagrangian, beta)
c4 = sym.diff(lagrangian, gamma)
mu_1 = m.addVar(name = "mu_1", lb = -gp.GRB.INFINITY, ub = gp.GRB.INFINITY)
mu_2 = m.addVar(name = "mu_2", lb = -gp.GRB.INFINITY, ub = gp.GRB.INFINITY)
lambda_1 = m.addVar(lb = 0, ub = gp.GRB.INFINITY)
lambda_2 = m.addVar(lb = 0, ub = gp.GRB.INFINITY)
beta = m.addVar(lb = -gp.GRB.INFINITY, ub = gp.GRB.INFINITY)
gamma = m.addVar(lb = -gp.GRB.INFINITY, ub = gp.GRB.INFINITY)
m.addConstr(c1==0)
m.addConstr(c2==0)
m.addConstr(c3==0)
m.addConstr(c4==0)
m.addConstr(f_lower <= mu_1*phi[z[0]] +mu_2*phi[z[1]])

m.setObjective(F_upper, gp.GRB.MINIMIZE)
m.optimize()
Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)
Thread count: 2 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 8 rows, 8 columns and 7 nonzeros
Model fingerprint: 0xe0107a48
Model has 2 quadratic objective terms
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e+00, 3e+00]
QObjective range [2e+00, 2e+00]
Bounds range [3e+00, 3e+00]
RHS range [1e+00, 2e+00]
Presolve time: 0.04s

Barrier solved model in 0 iterations and 0.04 seconds
Model is infeasible or unbounded
0

Comments

4 comments

Date Votes
  • Official comment
    Simranjit Kaur
    • Gurobi Staff
    Edited
    This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?.
  • Jaromił Najman
    • Gurobi Staff
    Edited

    Hi Varun,

    You are correct, Gurobi does not accept \(\texttt{sympy}\) objects as input for constraints. However, you can use the \(\texttt{sympy}\) objects to construct strings and work with these. For the sake of simplicity, let us assume that you want to add the partial derivatives of the term \(x_1 \cdot x_2 - x_1 + x_2\) to your model as \( \leq 0\) constraints, i.e., you want to add the constraints \( x_2 + 1 \leq 0\) and \( x_1 - 1 \leq 0\) to your model, which are computed by the \(\texttt{sympy}\) package.

    You can do so with the following code

    import gurobipy as gp
    from gurobipy import GRB
    import sympy as sym

    try:    
        # Create a new model
        m = gp.Model("test")

        # Create variable objects bound to the model m 
        m.addVar(vtype=GRB.CONTINUOUS,name="x1")
        m.addVar(vtype=GRB.CONTINUOUS,name="x2")
        
        # now the sympy part
        constraints = []
        variableNames = ["x1","x2"] # save all variable names
        x1 = sym.Symbol("x1")
        x2 = sym.Symbol("x2")

        test = x1*x2 - x1 + x2 # some test term
        # store the partial derivatives of the test term in the constraints list
        constraints.append(str(sym.diff(test,x1))) # = x2 + 1
        constraints.append(str(sym.diff(test,x2))) # = x1 - 1   	

        # write a Python file which is used to add the partial derivatives
        file = open("partialDerivatives.py","w")
        file.write("def addPartialDerivatives(model):\n")
        for v in variableNames:
        	file.write("    "+v+" = model.getVarByName(\""+v+"\")\n")
        for c in constraints:
            file.write("    model.addConstr("+c+" <= 0)\n")
        file.close()

        m.update() # essential to get all information
        # use the written file to add the partial derivative constraints
        import partialDerivatives
        partialDerivatives.addPartialDerivatives(m)

        # write an LP file to check whether everything worked
        m.write("myLP.lp")

    except gp.GurobiError as e:
        print('Error code ' + str(e.errno) + ': ' + str(e))

    except AttributeError:
        print('Encountered an attribute error')

    The above code writes a Python file \(\texttt{partialDerivatives.py}\) which holds the lines

    def addPartialDerivatives(model):
        x1 = model.getVarByName("x1")
        x2 = model.getVarByName("x2")
        model.addConstr(x2 - 1 <= 0)
        model.addConstr(x1 + 1 <= 0)

    In the file, we are first getting the correct variable objects via the getVarByName() function and then add the partial derivative constraints.

    This example is very simple and you will very likely have to adjust a few things for it to work with your mode. Please note, that the focus of the above code is not about performance or automation of each step but rather to present a possible solution.

    Best regards,
    Jaromił

    0
  • Satender .
    Edited
    • Gurobi-versary
    • First Question
    • First Comment

    Hi Jaromił Najman,

    I tried to extend your code above on matrix variables for adding KKT stationary condition but ran into an error. 

    import sympy as sym
    import numpy as np
    from sympy import Array
    from sympy import MatrixSymbol, Matrix
    from sympy import MatMul

    try:    
        # Create a new model
        m = gp.Model("test")

        # Example Data
        Mvar_rows, Mvar_cols = 1, 2
        c = np.array([1,10])
        A = np.array([
                      [1,2],
                      [3,-4]
        ])
        b = np.array([1,-2])

        # Variables objects for the model
        m.addMVar((Mvar_rows,Mvar_cols), vtype=GRB.CONTINUOUS,name="X")
        m.addMVar((1,A.shape[0]), vtype=GRB.CONTINUOUS, name = "L")

        # the sympy part
        X = MatrixSymbol('X', Mvar_rows, Mvar_cols)
        L = MatrixSymbol('L', 1, A.shape[0])         # Lagrange Multipliers

        constraints = []
        variableNames = ['X', 'L']
        Expression = c @ X.T
        for index in range(A.shape[0]):
          Expression += L[0,index]*(A[index,:]@X.T - b[index])

        # store the partial derivatives of the test term in the constraints list
        for i in range(Mvar_rows):
          for j in range(Mvar_cols):
            constraints.append(str(sym.diff(Expression[0], X[i,j]))) 

        # write a Python file which is used to add the partial derivatives
        file = open("partialDerivatives.py","w")
        file.write("def addPartialDerivatives(model):\n")
        for v in variableNames:
          file.write("    "+v+" = model.getVarByName(\""+v+"\")\n")
        for c in constraints:
            file.write("    model.addConstr("+c+" == 0)\n")
        file.close()
        m.update() # essential to get all information

        # use the written file to add the partial derivative constraints
        import partialDerivatives
        partialDerivatives.addPartialDerivatives(m)

        # write an LP file to check whether everything worked
        m.write("myLP.lp")

    except gp.GurobiError as e:
        print('Error code ' + str(e.errno) + ': ' + str(e))
    except AttributeError:
        print("Attribute Error")
    The file is written correctly (provided below).
    def addPartialDerivatives(model):
        X = model.getVarByName("X")
        L = model.getVarByName("L")
        model.addConstr(L[0, 0] + 3*L[0, 1] + 1 == 0)
        model.addConstr(2*L[0, 0] - 4*L[0, 1] + 10 == 0)

     The error is with getVarByName used on matrix variables. Error is given below, the program is reading var names L and X as none types.

    TypeError                                 Traceback (most recent call last)


    <ipython-input-155-2c179d7a7531> in <module>()
    
     49     # use the written file to add the partial derivative constraints
    
         50     import partialDerivatives
    
---> 51     partialDerivatives.addPartialDerivatives(m)
     
         52 
    
     53     # write an LP file to check whether everything worked



    /content/partialDerivatives.py in addPartialDerivatives(model)
      
          2     X = model.getVarByName("X")
    
      3     L = model.getVarByName("L")
    
----> 4     model.addConstr(L[0, 0] + 3*L[0, 1] + 1 == 0)
    
      5     model.addConstr(2*L[0, 0] - 4*L[0, 1] + 10 == 0)




    TypeError: unsupported operand type(s) for -: 'NoneType' and 'int'

     

    0
  • Jaromił Najman
    • Gurobi Staff
    Edited

    Hi,

    Your variables are named \(\texttt{L[0,i]}\). Thus, you have to use these exact names to retrieve them when using the getVarByName method. Something like the following should work

    L = []
    for i in A.shape[0]:
      L.append(model.getVarByName("L[0,%d]"%i))

    Note that you might have to pass the \(\texttt{A.shape[0]}\) value to your function.

    Best regards, 
    Jaromił

    1

Post is closed for comments.

Didn't find what you were looking for?

New post