Solving Non-linear optimization in Gurobi
AnsweredHello, I am a new Gurobi user. I am wondering if it is possible to solve a problem like in this photo in the Gurobi student version.
-
-
Gurobi Staff
Hi.
It is possible as long as \(\epsilon_{\beta} > 0\). You will have to introduce multiple auxiliary variables. First you have to formulate the term \((lp + \epsilon_{\alpha} \cdot \sum X_{ir})^{\epsilon_{\beta}}\) as
\[\begin{align*}
z_{1i} &= lp + \epsilon_{\alpha} \cdot \sum X_{ir}\\
z_{2i} &= z_1^{\epsilon_{\beta}}
\end{align*}\]The power term can be implemented with the addGenConstrPow method.
You can then reformulate the multilinear term \(VOT_i \cdot t_f \cdot z_{2i} \cdot X_{ir}\) as described in How do I model multilinear terms in Gurobi?
Please note that since the expression you want to optimize is highly nonlinear, it is necessary that you provide tight bounds for every participating variable in order to avoid numerical issues.
Best regards,
Jaromił1 -
-
Hi Jaromił,
All the variables are binary. I am not sure how to put tight bounds on binary variables. Could you please provide a simple example?
0 -
-
-
All the variables are binary. I am not sure how to put tight bounds on binary variables. Could you please provide a simple example?
If all variables are binary, then you should not have to provide any additional bound information.
1 -
-
Hi Jaromił,
I reformulated the problem as you suggested using addGenConstrPow method. But I kept the multilinear terms as they were. The model resulted in 0 objective value solution which is wrong. I am wondering if this is because I did not convert the multilinear terms.
Also, in the objective function, since all the terms except \(X_{ir}\) are given as parameter values (not decision variables) do I really need to convert the multilinear terms as I suggested? I read the article "How do I model multilinear terms in Gurobi" but I am not sure how would I do this in my case. I am confused because each term in my objective function has different domains. So while I do understand that I can set \(x*y = w\) but I am not sure about how would I go to model \(x_i \times y_j = w, \ i \ne j\). What would be the dimension of \(w\)?
A reproducible version of my code is given below. Ideally, I want the objective function in my code. But I am not sure how I would convert the multilinear terms.
For some reason when I paste the code in the code-box, all codes cram into a single line. I put the code in the link below as well:
https://drive.google.com/file/d/1U7J-QpZcG5kKKUOGdgA8VSDXLXYlnrlG/view?usp=sharing
import gurobipy as gp
from gurobipy import Model, GRB, quicksum
customer_VOT_set = {'2_OD1': 2.1, '10_OD2': 1.8, '1_OD1': 2.1, '3_OD2': 1.8, '3_OD1': 2.1, '5_OD2': 1.8}
pot_route = {'2_OD1': [124, 1234, 134], '10_OD2': [24, 234], '1_OD1': [24, 234], '3_OD2': [24, 234], '3_OD1': [124, 1234, 134], '5_OD2': [24, 234]}
link_in_route = {34: [134, 234, 1234], 12: [1234, 124], 13: [134], 23: [234, 1234], 24: [24, 124]}
which_links_in_route = {134: [13, 34], 234: [23, 34], 1234: [12, 23, 34], 24: [24], 124: [12, 24]}
road_travel_time = {12: 4.5, 13: 5, 24: 6, 23: 0.8, 34: 5}
all_unique_routes = [134, 234, 1234, 24, 124]
all_unique_links = [34, 12, 13, 23, 24]
alpha = 0.15
pass_vehicles = 2200
capacity = 2400
mdl = Model('OTR')
X = {}
for item in pot_route:
for element in pot_route[item]:
X[item,element] = mdl.addVar(vtype=GRB.BINARY, name="X[{0},{1}]".format(item,element))
mdl.addConstrs((quicksum(X[item,j] for j in pot_route[item])==1) for item in pot_route)
L = mdl.addVars((item for item in road_travel_time), name = "L")
ZC_rep = {}
ZC = {}
LC = {}
for item in all_unique_links:
ZC_rep[item] = mdl.addVar(lb=0, name="ZC_rep[{0}]".format(item))
route_set = []
for element in which_links_in_route:
if item in which_links_in_route[element]:
route_set.append(element)
mdl.addConstr(ZC_rep[item]==(pass_vehicles+3*(quicksum(X[i,j] for i in pot_route for j in route_set if j in pot_route[i])))/capacity, name="Z_rep[{0}]".format(item))
for item in all_unique_links:
LC[item] = mdl.addVar(lb=-GRB.INFINITY, name="LC[{0}]".format(item))
ZC[item] = mdl.addVar(lb=0, name="ZC[{0}]".format(item))
mdl.addGenConstrPow(ZC_rep[item],ZC[item],4,"rep","FuncPieces=1000")
for item in all_unique_links:
mdl.addConstr(LC[item]==road_travel_time[item]*(1+alpha*ZC[item]), name="Link_cost[{0}]".format(item))
RC = {}
for item in all_unique_routes:
RC[item] = mdl.addVar(lb=0, name="RC[{0}]".format(item))
mdl.addConstr(RC[item]== quicksum(L[i] for i in which_links_in_route[item]), name="RC_cost[{0}]".format(item))
mdl.write("myLP.lp")
mdl.setObjective(quicksum(customer_VOT_set[item]*X[item,element]*RC[element]
for item in customer_VOT_set for element in all_unique_routes if element in pot_route[item]),GRB.MINIMIZE)
#mdl.params.NonConvex = 2
mdl.update()
mdl.optimize()0 -
-
-
Hi Tanvir,
Also, in the objective function, since all the terms except \(X_{ir}\) are given as parameter values (not decision variables) do I really need to convert the multilinear terms as I suggested?
No, in this case you don't have to reformulate these terms. It is only necessary if you multiply more than 2 optimization variables with each other.
The model resulted in 0 objective value solution which is wrong. I am wondering if this is because I did not convert the multilinear terms.
Do you know what the optimal objective value should be? Do you know the solution point values?
Best regards,
Jaromił0 -
-
Hi Jaromił,
I only know that due to the constraints, any feasible solution cannot have 0 objective value. I can also figure out a feasible solution for this small data but not for my actual one which is large.
0 -
-
-
Hi Tanvir,
Not having a solution at hand makes it very hard to help. You could analyze the LP file you write to make sure that all constraint look like what they should. Please also check all variable bounds. Note that if a variable is not listed as free (or with explicit bounds) then it has a lower bound of 0 and an infinite upper bound. Note that you should write the LP after your set your objective. Please also note that the update call is not needed.
mdl.setObjective(quicksum(customer_VOT_set[item]*X[item,element]*RC[element]
for item in customer_VOT_set for element in all_unique_routes if element in pot_route[item]),GRB.MINIMIZE)
#mdl.params.NonConvex = 2
mdl.write("myLP.lp")
mdl.optimize()You can also write the optimal solution found by Gurobi to a file
mdl.optimize()
if mdl.SolCount > 0:
mdl.write("solution.sol")and inspect the solution point values. Maybe this helps.
Best regards,
Jaromił1 -
-
Hello Jaromił,
Your suggestions have been helpful and I finally found a solution. Thank you.
Just out of curiosity, what algorithm Gurobi uses to solve integer (in my case non-linear) programming problems?
0 -
-
-
Hi Tanvir,
Just out of curiosity, what algorithm Gurobi uses to solve integer (in my case non-linear) programming problems?
Pure integer problems can be reformulated to MIPs via linearization techniques, see our webinar by Ed Klotz on Models with products of binaries. If your model also has products with continuous variables, then Gurobi uses a spatial B&B approach. You can find more information in our nonconvex webinar and nonconvex tech talk.
Best regards,
Jaromił0 -
Please sign in to leave a comment.
Comments
9 comments