Is there any api which can be used to retrieve the standard form of a linear program in python.
AnsweredWhile solving large-scale optimization problem, I want to know if Gurobi provides a API which can obtain the coefficient matrix (A,b,C,Aineq,bineq ) of a LP problem.
min C'X
s.t. Ax=b
Aineqx<=bineq
-
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,
Matthias1 -
Hi Matthias
There is still a problem, can I retrieve A & Aineq respectively.
Best
Canqi
0 -
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
Matthias0 -
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 -
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 -
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
Matthias0 -
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 -
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,
Matthias0
Please sign in to leave a comment.
Comments
8 comments