• 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

Hi Matthias

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

Best

Canqi

• 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

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

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 VRPimport mathimport randomfrom gurobipy import *from scipy.io import loadmat, savematimport numpy as npimport sysimport csvimport timeclass barX:    passclass Psi_fc:    passclass Psi_ic:    pass# Initial paras# map = 100# node=32def 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,xif __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')
• 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

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

• 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