Solving Nonlinear Integer model with Gurobi
OngoingHello,
I am trying to implement a single optimization vehicle platooning problem in the most straightforward way to approach a nonlinear integer program (INLP).
The result I got from the model does not make any sense. I expect my solution after some iteration.
Could you check what is the problem with my model? the model source from paper KTM page 13
fuel_saving = 0.1
nodes=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
time = [12, 24, 20, 13, 11, 18, 19, 28, 15, 21] #the earliest pickup time at node
edges = [(0, 8),(0, 4),(0, 7),(0, 1),(0, 3),(1, 8),(1, 7),(1, 5),(1, 2),(2, 8),(2, 5),(2, 7),(2, 3),(3, 7),(3, 5),
(3, 8),(3, 4),(4, 5),(4, 9),(5, 7),(7, 9),(7, 8)]
distance = {(0, 8): 5,(0, 4): 3,(0, 7): 1,(0, 1): 10,(0, 3): 7,(1, 8): 2,(1, 7): 8,
(1, 5): 3,(1, 2): 3,(2, 8): 6,(2, 5): 10,(2, 7): 7,(2, 3): 4,(3, 7): 3,
(3, 5): 8,(3, 8): 5,(3, 4): 7,(4, 5): 7,(4, 9): 3,(5, 7): 6,(7, 9): 5,(7, 8): 7}
trucks = [0,1]
vehicles_destination=[5, 1]
vehicles_origin=[1, 8]
mdl = Model('CVPP')
x = {}
for truck in trucks:
x[truck] = mdl.addVars(time,edges,vtype=GRB.BINARY)
delta= mdl.addVars(time, edges,vtype=GRB.BINARY, name="delta")
Minimization of cost
\[\begin{align} Min \ Z \ = \sum_{k}x^k_{t,ij}  \delta_{t,ij}\eta \left ( \sum_{k} x^k_{t,ij} 1 \right ) \ \forall i,j \ \epsilon G; e_{i,j} \ \epsilon E ; t \ \epsilon T \end{align}\]
obj1 = LinExpr()
f1 = mdl.addVar()
for truck in trucks:
obj1=quicksum(x[truck][t,i,j] for i, j in edges for t in time)
obj2 =(delta.sum("*","*","*")*fuel_saving)*(obj1 1)
c_f1 =mdl.addConstr(obj1obj2 == f1)
O_f1 = mdl.setObjective(f1,GRB.MINIMIZE)
for truck in trucks:
mdl.addConstr(delta.sum("*","*","*")<= quicksum(x[truck][t,i,j] for i, j in edges for t in time))
for truck in trucks:
for t in range(len(time)  1):
con_1 = mdl.addConstr(quicksum(x[truck][time[t],i,j] for i,j in edges) == quicksum(x[truck][time[t+1],i,j] for i,j in edges))
for truck in trucks:
for t in time:
for i,j in edges:
if vehicles_origin[k] == i:
rhs = x[truck][t,i,j]
con_2 = mdl.addConstr(quicksum(x[truck][t,i,j] for i,j in edges if t == time[1] ) == rhs)
Please help with this model

Hi John,
1. The decision variables are defined correctly. Alternatively, you can define them as below to remove the outer \(\texttt{for}\) loop. To represent the set \(T\) in the paper which is the set of equidistant time periods, let us use \(\texttt{time_steps}\) instead of \(\texttt{time}\) where \(\texttt{time_steps}\) = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9].
# x_k_t_ij = 1 if truck k travels between time steps t and t+1 along the edge (i, j)
x = model.addVars(trucks, time_steps, edges, vtype=GRB.BINARY, name="x")
# delta_ij
delta = model.addVars(time_steps, edges, vtype=GRB.BINARY, name="delta")2. The objective function, shown below, is not modelled correctly and needs to be rewritten completely.
\[\sum_t \sum_{ij \in e} c_{t, ij} = \sum_t \sum_{ij \in e} f_{ij} \big(\sum_k x^k_{t, ij}  \eta \delta_{t, ij} (\sum_k x^k_{t, ij} 1)\big)\]
# Objective function
objective = gurobipy.quicksum(
distance[(i, j)]
* (
x.sum("*", t, i, j)
 fuel_saving * delta[t, i, j] * (x.sum("*", t, i, j)  1)
)
for t in time_steps
for (i, j) in edges
)
model.setObjective(objective, sense=GRB.MINIMIZE)3. Constraints need to be rewritten as well.# Constraint 3.1
model.addConstrs(
(
x.sum(k, t, "*", i) == x.sum(k, t + 1, i, "*")
for k in trucks
for t in time_steps[:1]
for (i, _) in edges
),
name="constr_3_1",
)
# Constraint 3.2
model.addConstrs(
(x.sum(k, 0, vehicles_origin[k], "*") == 1 for k in trucks),
name="constr_3_2",
)
# Constraint 3.3
T = len(time_steps)  1
model.addConstrs(
(x.sum(k, T, "*", vehicles_destination[k]) == 1 for k in trucks),
name="constr_3_3",
)
# Constraint 3.5
model.addConstrs(
(delta[t, i, j] <= x.sum("*", t, i, j) for t in time_steps for (i, j) in edges),
name=constr_3_5"
)Please be aware that given the current input, the model is infeasible. The origin node for truck 2 is 8 and the destination node is 1. However, in the set of edges, there is no edge leaving node 8 and there is no possible route for truck 2 to travel from node 8 to node 1.
Best regards,
Maliheh
1 
Hello Maliheh,
Thank you for your answer is very helpful.
In my logic, my graph(input) is bidirectional edges. that means there is a route from 8 to 1.
But I change input edges and distance like this.edges = [(i, j) for i in nodes for j in nodes if i != j]
distance = {(i, j): random.randint(5, 10) for i, j in edges}time input I add to this model because I need to a constrain that ensure each truck reach their destination before the deadline.
Now I am so confused because of time_steps input.t_max = [80,120] #deadline of each truck
so now I write code like that remove time input and add t_ij parameter
time_ij = {}
for i,j in edges:
time_ij[i,j] = distance[i,j]/v
for truck in trucks:
mdl.addConstr(x.sum("*", t, i, j)*time_ij[i,j]for i,j in edges) <= t_max[truck])
Could you please tell me how can I add this constrain?Thank you
0 
Hi John,
Sure, the graph can be bidirectional. We just need to make sure that the model knows about it :) What you did sounds alright to me.
Defining \(t_{ij}\) as the travelling time from node \(i\) to node \(j\) and \(t^k_{\mathrm{max}}\) as the total possible travelling time for truck \(k\), I assume the new constraint is:
\[\sum_{ij} t_{ij} \sum_{t} x^k_{t, ij} \leq t^k_{\mathrm{max}}, ~~ \forall k \].
You can then implement it as:
model.addConstrs(
(
gurobipy.quicksum(x.sum(k, "*", i, j) * time[(i, j)] for (i, j) in edges)
<= t_max[k]
for k in trucks
),
name="constr_travelling_time",
)
Best regards,Maliheh
1 
Hello,
Thank you for your help.constrain for travelling time doesn't make any changes to my model. something I miss in my model
for example, if track1 max time is 20 and start node 1 and end node is 5.
possible path 1: {(1,2):5, (2,3):5, (3,4):6, (4,5):5}
2:{(1,2):5, (2,3):5, (3,5):8}
path 1 total time is 21
path 2 total time is 18
select path 2 using time constrain:
moreover, I change for getting complete graph
edges = [(i, j) for i in nodes for j in nodes if j>i]
after changing to complete graph model doesn't have any solution.
Thank you.0 
Hi John,
One idea to debug your model is to write it into an LP format using the Model.write() method. You can then open the file using any text editor to figure out what it is missing in the model.
Best regards,
Maliheh
0 
Hi Maliheh,
After so much reading and a lot of tries, I thought of covert this constrain to the objective function.Currently, I have 2 objectives for my model I need to use the epsilon constrain method to optimize the solution
use networkx package for creating graph
minimize total time of all truck in system\[\begin{align} min = \sum_h x^h_{i,j} \sum_e t_{i,j} \end{align}\]
h = truck
objective_function_main_second = quicksum(x[edge][h] * quicksum(time[edge] for edge in road_network.edges())for h, locations in group.items())
model.setObjective(objective_function_main_second, GRB.MINIMIZE)
I modified variables x and yx = {}
y = {}
for edge in road_network.edges():
x[edge] = {}
y[edge] = model.addVar(vtype=GRB.BINARY)
for h, locations in group.items():
x[edge][h] = model.addVar(vtype=GRB.BINARY)
Could you please check my objective function is correct?0
Please sign in to leave a comment.
Comments
6 comments