• Gurobi Staff

Hi Tanay,

One way you can debug this is to use an interactive Python session (such as a Jupyter notebook) and set a value for k, and see what results:

k = 0np.sum(x[k],axis(0)) == X0[k]

Then you can drill down, inspecting smaller components, such as np.sum(x[k]), and see what they result in.

Note the axis(0) is not valid, as there is no object axis defined. Perhaps you meant

np.sum(x[k], axis=0)

or

np.sum(x[k], 0)

However note that in general you can't pass Gurobi objects (like the MVar x[k]) to numpy functions - numpy does not know what to do with these things.  Instead we have a number of functions defined on the MVar class, like MVar.sum(), as well as some defined in the gurobipy namespace: Matrix-Friendly API Functions.

I expect you will want to use the following:

x[k].sum(axis=0) == X0[k]

Lastly, it looks like you are trying to add multiple constraints at once (for k in range(3)), in which case you need to use Model.addConstrs not Model.addConstr, or alternatively loop:

for k in range(3):     model.addConstr(x[k].sum(axis=0) == X0[k])

- Riley

Hello Riley and Gurobi team,

Thank you for your suggestions. I implemented a version of the Miller-Tucker-Zemlin constraints because in the multiple classes of vehicles, the solutions had subtours. However, now the solutions I'm getting are in clear violations of the constraints.

I'm getting a solution like x[0][0][i] = 1 for all i belonging to (1,2,3...7)

By the constraint that if j!=i (lets keep j as 1 and i as 0), then

demTillNode for the 1st node should be >= demTillNode for [0] + x[k][0][1] * dem of 1st node (which is 1.2)

So demTillNode for 1st node should've been >= 1.2 at least, but here you can see that it's 0. Could someone help me in what's going wrong? Thank you!

!pip install gurobipyimport gurobipy as gpparams = {"WLSACCESSID":"xxxxxxxxxxxxx","WLSSECRET":"xxxxxxxxxxxx","LICENSEID":xxxxx,}env = gp.Env(params=params)model= gp.Model(env=env)from gurobipy import GRBimport math#K is the count of Vehicle TypesK= 3#N is the universe of Nodes with [0] being the DCN = 8#s is the matrix of the distances between any two nodes on the graph. First row & column (index 0) are measured from the DCs = [[0.0,60.12,62.48,64.81,106.91,74.37,86.4,69.21],[60.12,0.0,13.91,4.76,52.14,15.62,29.3,9.34],[62.48,13.91,0.0,15.25,60.72,24.75,37.97,18.19],[64.81,4.76,15.25,0.0,47.87,11.14,24.93,4.6],[106.91,52.14,60.72,47.87,0.0,36.77,22.98,43.56],[74.37,15.62,24.75,11.14,36.77,0.0,13.8,6.82],[86.4,29.3,37.97,24.93,22.98,13.8,0.0,20.57],[69.21,9.34,18.19,4.6,43.56,6.82,20.57,0.0]]#DMax is the carrying capacities of the 3 vehicle types in cbmDMax = [3.0,6.0,15.0]#c is the matrix of converting each distance to currency (e.g.- PHP) based on the vehicle type deployed between those two nodes.c = [10.0,15.0,35.0]#cnew is the cost of deploying any new vehiclecnew = 2000.0#d is the demand per node of the graphd = [0,1.2,1.5,4.9,0.5,1.7,1.2,1.0]#M is the max number of nodes allowed on one ride of one vehicleM = 2.0#L is the max distance which any vehicle can cover in a single ride (i.e. in one day) due to urban traffic congestionL = 150# x[k][i][j] defines whether a vehicle of the kth type traversed from node i to node jx = [[[ [] for _ in range(N)] for _ in range(N)] for _ in range(K)]for k in range(K):for i inrange(N):for j inrange(N):x[k][i][j] = model.addVar(vtype=GRB.BINARY, name= f'x[{k}][{i}][{j}]')#Need z to multiply y with cost per distance as per each vehicle typez = [[[ [] for _ in range(N)] for _ in range(N)] for _ in range(K)]for k in range(K):for i inrange(N):for j inrange(N):z[k][i][j] = model.addVar(name= f'z[{k}][{i}][{j}]')lenTillNode = [[[] for _ in range(N)] for _ in range(K)]distTillNode = [[[] for _ in range(N)] for _ in range(K)]demTillNode = [[[] for _ in range(N)] for _ in range(K)]for k in range(K):for i inrange(N):lenTillNode[k][i] = model.addVar(name=f'lenTillNode[{k}][{i}]')distTillNode[k][i] = model.addVar(name=f'distTillNode[{k}][{i}]')demTillNode[k][i] = model.addVar(name=f'demTillNode[{k}][{i}]')#No ride from the same node to same node and no ride to the DCfor k in range(K):for i inrange(N):for j inrange(N):model.addConstr(x[k][i][i] == 0, name='Intranode travel set as zero')model.addConstr(x[k][i][0] == 0, name='The DC has no receiving')model.addConstr(x[k][i][j] + x[k][j][i] <= 1, name='No backhaul or looping allowed')#For all nodes which are not 0 (i.e. the DC), the total indexes incoming are 1, so that 1 vehicle of a certain type from 1 node visits this nodefor j in range(1,N):model.addConstr(sum(x[k][i][j] for i in range(N) for k in range(K)) == 1, name='Each node touched just once')#Define vehicle-type wise variable costs based on distance travelledfor k in range(K):for i inrange(N):for j inrange(N):model.addConstr(x[k][i][j] * s[i][j] * c[k] == z[k][i][j], name='Binary edges assigned with costs')#For any edge, a vehicle of type-k can traverse iff a vehicle of the same type has reached. We assume that vehicles can't change their type at nodes.#for k in range(K):#for i in range(1,N):#model.addConstr(sum(x[k][i][j] for j in range (1,N)) <= sum(x[k][jdash0][i] for jdash0 in range(N)),name='Conservation of vehicles at node')#model.addConstr(sum(x[k][i][j] for j in range(1,N)) <= 1, name='Max one vehicle exits node')for k in range(K):model.addConstr(lenTillNode[k][0] == 0)model.addConstr(distTillNode[k][0] == 0)model.addConstr(demTillNode[k][0] == 0)for i inrange(1,N):model.addConstr(lenTillNode[k][i] >= 1)model.addConstr(lenTillNode[k][i] <= M)model.addConstr(distTillNode[k][i] <= L)model.addConstr(demTillNode[k][i] <= DMax[k])for j inrange(1,N):if j != i:model.addConstr(lenTillNode[k][j] >= lenTillNode[k][i] + x[k][i][j])model.addConstr(distTillNode[k][j] >= distTillNode[k][i] + x[k][i][j] * s[i][j])model.addConstr(demTillNode[k][j] >= demTillNode[k][i] + x[k][i][j] * d[j])model.setObjective(sum(z[k][i][j] for i in range(N) for j in range(N) for k in range(K)) + cnew * sum(x[k][0][i] for i in range(N)), GRB.MINIMIZE)model.optimize()sysNames = model.getAttr('VarName',model.getVars())sysValues = model.getAttr('X',model.getVars())for idx in range(len(sysNames)):print(sysNames[idx],' ',sysValues[idx])

• Gurobi Staff

Hi Tanay,

My advice is to write your model to a LP file, and open the file in a text editor.  You will want to verify that the constraint you believe to be violated is present in the model.  You can then manually calculate the values of the LHS and RHS using the solution values and verify the constraint is violated.  I expect one of these two steps will fail, but the learning will help fix the model.

- Riley