Skip to main content

Is there any api which can be used to retrieve the standard form of a linear program in python.

Answered

Comments

8 comments

  • Matthias Miltenberger
    Gurobi Staff Gurobi Staff

    Hi Canqi,

    You can use the Python method model.getA() to get the entire constraint matrix in one call. The objective coefficients and the right-hand sides can be retrieved like this:

    m.getAttr('Obj', m.getVars())
    m.getAttr('RHS', m.getConstrs())

    The constraints' senses can be retrieved accordingly:

    m.getAttr('Sense', m.getConstrs())

    Cheers,
    Matthias

    1
  • Canqi Yao
    Gurobi-versary
    Conversationalist
    Curious

    Hi Matthias 

    There is still a problem, can I retrieve A & Aineq respectively.

    Best

    Canqi

    0
  • Matthias Miltenberger
    Gurobi Staff Gurobi Staff

    Hi Canqi,

    You can use the constraints' senses to distinguish inequalities and equalities. You may also go through the constraints one by one and build up \(\tt{A}\) and \(\tt{Aineq}\) iteratively based on the sense.

    There is no dedicated method to get only the equations' coefficients nor only the inequalities' coefficients.

    Cheers
    Matthias

    0
  • Canqi Yao
    Gurobi-versary
    Conversationalist
    Curious

    Hi Matthias 

    Using your recommended code cannot obtain the correct right hand side value for inequality constraints, which has been verified by hand. Can you explain the mechanism of adding constraints with Gurobi. 

    m.getAttr('RHS', m.getConstrs())

    Best

    Canqi

     

    0
  • Canqi Yao
    Gurobi-versary
    Conversationalist
    Curious

    As shown in the following figure, the 1st element of bineq should be 5000 in this case ( T[i,i]==0 for i in range(N) ). However, running the code gives me a wrong RHS value.

     

     

    # branch and cut used in VRP
    import math
    import random
    from gurobipy import *
    from scipy.io import loadmat, savemat
    import numpy as np
    import sys
    import csv
    import time


    class barX:
    pass


    class Psi_fc:
    pass


    class Psi_ic:
    pass


    # Initial paras

    # map = 100
    # node=32
    def BVRP(vehicle, node):
    ###### Main function
    # node=32
    # vehicle=3
    data = loadmat('{}vehicles/rc101_{}_test.mat'.format(vehicle, node))
    x_sol = loadmat("x_sol.mat")
    x = x_sol['x_sol']
    # print(data)



    Emax = data['Emax'][0][0]
    rmax = Emax
    e = data['e']
    T = data['T']
    M = int(data['M'])
    p = data['p'][0]
    g = data['g'][0]
    t = data['t'][0]
    c = data['c'][0]
    K = int(vehicle)
    N = int(data['N'])
    soc = 0.2

    ite = 1
    tol = 1e-1
    psi_ic = []
    psi_fc = []
    UB = []
    LB = []
    UB.append(1e7)
    LB.append(-1e7)
    s = time.time()
    print('x[0,0,0]', x[0,1,0],t[0],T[0,0])

    master_method = 'fast'
    "Initial solution"
    m1 = Model()
    m1.setParam('InfUnbdInfo', 1)
    E1 = m1.addVars(N, K, vtype=GRB.CONTINUOUS, name="E1")
    E2 = m1.addVars(N, K, vtype=GRB.CONTINUOUS, name="E2")
    eta = m1.addVars(N, N, K, vtype=GRB.CONTINUOUS, name="eta") # ****************
    epsilon = m1.addVars(N, N, K, vtype=GRB.CONTINUOUS, name="epsilon") # ****************
    r = m1.addVars(N, K, vtype=GRB.CONTINUOUS, name="r")

    m1.addConstrs(( 0 >= -t[j] + T[i, j] + g[i] * r[i, k] + t[i] - M * (1 - x[i, j, k]) for i in range(N) for j in range(N) for k in range(K)),name='c')
    # print(m1.constr)
    m1.addConstrs(
    0 <= M * (1 - x[i, j, k]) + E1[i, k] - e[i, j] * x[i, j, k] + r[i, k] - E2[j, k] for i in range(N) for j in
    range(N) for k in range(K))
    m1.addConstrs(
    E1[i, k] - e[i, j] * x[i, j, k] + r[i, k] - E2[j, k] - M * (1 - x[i, j, k]) <= 0 for i in range(N) for j in
    range(N) for k in range(K))

    m1.addConstrs(
    eta[i, j, k] >= g[i] * r[i, k] - M * (1 - x[i, j, k]) for i in range(N) for j in range(N) for k in range(K))
    m1.addConstrs(
    epsilon[i, j, k] >= p[i] * r[i, k] - M * (1 - x[i, j, k]) for i in range(N) for j in range(N) for k in range(K))
    ## soc and bilinear
    m1.addConstrs(E1[j, k] <= Emax for j in range(N) for k in range(K))
    m1.addConstrs(r[j, k] <= rmax for j in range(N) for k in range(K))
    m1.addConstrs(E2[j, k] <= Emax for j in range(N) for k in range(K))
    m1.addConstrs(eta[i, j, k] >= 0 for i in range(N) for j in range(N) for k in range(K))
    m1.addConstrs(epsilon[i, j, k] >= 0 for i in range(N) for j in range(N) for k in range(K))
    m1.addConstrs(E1[j, k] >= 0 for j in range(N) for k in range(K))
    m1.addConstrs(E2[j, k] >= 0 for j in range(N) for k in range(K))
    m1.addConstrs(r[j, k] >= 0 for j in range(N) for k in range(K))
    m1.addConstrs(E1[0, k] == soc * Emax for k in range(K))
    m1.addConstrs(E1[j, k] == E2[j, k] for j in range(N) for k in range(K))

    # m1.addConstrs(x[i,j,k]<=1 for k in range(K) for i in range(N) for j in range(N))
    # m1.addConstrs(0<=x[i,j,k] for k in range(K) for i in range(N) for j in range(N))
    # m1.addConstrs(x.sum( i,'*', '*') <= 1 for i in range(N) if i != 0)
    #
    # m1.addConstrs(x[i, i,k] == 0 for i in range(N) for k in range(K))
    # m1.addConstrs(x[ N - 1, i,k] == 0 for i in range(N) for k in range(K))
    # m1.addConstrs(x[ i, 0,k] == 0 for i in range(N) for k in range(K))
    #
    # # Flow conservation
    # m1.addConstrs(x.sum( 0, '*',k) - x.sum( '*', 0,k) == 1 for k in range(K))
    # m1.addConstrs(x.sum( N - 1, '*',k) - x.sum('*', N - 1,k) == -1 for k in range(K))
    # m1.addConstrs(x.sum( i, '*',k) - x.sum('*', i,k) == 0 for i in range(N - 1) if i != 0 for k in range(K))

    m1.setObjective(quicksum(eta) + quicksum(epsilon)) # +epsilon
    e_time = time.time()
    m1.optimize()

    A = m1.getA()
    # f=m.getCoeff()
    f = m1.getAttr('Obj', m1.getVars())
    b = m1.getAttr('RHS', m1.getConstrs())
    sense = m1.getAttr('Sense', m1.getConstrs())
    # count element '='
    count= sense.count('=')
    counine=sense.count('<')
    print(len(sense) - count)
    Aineq = A[0:len(sense) - count, :]
    Aeq = A[len(sense) - count:len(sense), :]
    bineq = b[0:len(sense) - count]
    beq = b[len(sense) - count:len(sense)]
    # np.transpose

    if m1.status == GRB.status.OPTIMAL:
    print('Optimal solution')
    obj = m1.getObjective()
    obj_sol = obj.getValue()
    # dualbd=
    else:
    print('Unbounded')
    obj_sol= math.inf
    dualbd = m1.getAttr('FarkasDual', m1.getConstrs())

    # m1.getAttr('UnbdRay')


    # print( "(A) Matrix coeff. %s" %Aeq)
    # print( "(A) Matrix coeff. %s" %f)
    # print( "(A) Matrix coeff. %s" %b)
    # print( "(A) Matrix coeff. %s" %sense)
    return Aeq, Aineq, beq, bineq, f, sense,obj_sol,m1,A,b,x

    if __name__ == '__main__':
    [Aeq, Aineq, beq, bineq, f, sense,obj_sol,m1,A,b,x] = BVRP(vehicle=3, node=31)
    print(min(beq))
    savemat('coeff.mat', mdict={'Aeq': Aeq, 'A': Aineq, 'beq': beq, 'b': bineq, 'f': f})
    # savemat('dualbd.mat',mdict={'dualbd':dualbd,'obj_sol':obj_sol,'m1':m1})

    print('Saved data')
    0
  • Matthias Miltenberger
    Gurobi Staff Gurobi Staff

    Hi Canqi,

    There are constants in the constraint definition:

    m1.addConstrs((0 >= -t[j] + T[i, j] + g[i] * r[i, k] + t[i] - M * (1 - x[i, j, k])
    for i in range(N) for j in range(N) for k in range(K)), name='c')

    The only variables in this constraint definition are \(r_{i,k}\) - everything else are just constant terms that add up to the value you are seeing in the output.

    Please note that constant values are always moved to the right-hand side to simplify the constraint.

    Cheers
    Matthias

    0
  • Canqi Yao
    Gurobi-versary
    Conversationalist
    Curious

    Hi Matthias 

    Note that when i=0, j=0, and k=0, this constraint =

    0 >= -t[0]+ T[0,0] + g[0] * r[0,0] + t[0] - M * (1 - x[0,0,0]

    which can be further simplified for T[0,0]=0

    -g[0]*r[0,0]>= - M * (1 - x[0,0,0]

    Due to x[0,0,0]=0, the correct right hand side value of the above constraint should be -M (M=5000).

    Could you please explain this problem?

    Best 

    Canqi

    0
  • Matthias Miltenberger
    Gurobi Staff Gurobi Staff

    Hi Canqi!

    There is no way for me to check your claim. Gurobi should not modify the data in such a way. I suggest you write down a very simple model that only contains the first constraint and then you can better check what's happening.

    Another approach is to write out the LP file using \(\texttt{write("model.LP")}\) after the model construction is finished. Then you can check whether the generated model matches your mathematical model.

    Cheers,
    Matthias

    0

Please sign in to leave a comment.