trajectory planniing problem
AnsweredI am working on the trajectory planning with a goal minimizing the fuel consumption whose function is non linear. I have crafted the following codes according to the model descriptions https://1drv.ms/w/s!AuGOQ6JTPZhskVk5C70tBo9C0Xtm?e=uEMJsu . The following codes turns out infeasible. Can anyone help me on this one? Thank you for any help.
def OptimalTrajectory1(timesteps=18, initialSpeed=25, distance=300, minimumSpeed=1, maximumSpeed=30):
### Zhihong Yao integrated ....
# Parameters (constants)
N = timesteps # Instants or Steps
delT = 1 # resolution
r1 = 1.0
u_max_x = 1
# state and input variables initialization (variables)
x0 = 1.0
vx0 = initialSpeed
ux0 = 0.1
# final state to reach at:
xN = distance
vxN = 18 # final speed when the vehicle arrive at the intersection
solutions = {}
alpha = 0.666
beta1 = 0.072
beta2 = 0.0344
d1 = 0.269
d2 = 0.0171
d3 = 0.000672
mass = 1680
bigM = 5000
# Create a new model
m = Model("lqrmodel")
m.Params.OutputFlag = 0
# Create variables
x = m.addVars(N, lb=0.0, ub=float('inf'), vtype=GRB.CONTINUOUS, name="x")
vx = m.addVars(N, lb=minimumSpeed, ub=maximumSpeed, vtype=GRB.CONTINUOUS, name="vx")
ux = m.addVars(N, lb=-u_max_x, ub=u_max_x, vtype=GRB.CONTINUOUS, name="ux")
fuelBinary = m.addVars(N, vtype = GRB.BINARY, name="FCIndicator")
power = m.addVars(N, lb=float('inf'), ub=float("inf"), vtype=GRB.CONTINUOUS, name="power")
threeOrder = m.addVars(N, lb = 0.001, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
twoOrder = m.addVars(N, lb = 0.001, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
fuelConsumption = m.addVars(N, lb= 0.00, ub = float('inf'), vtype = GRB.CONTINUOUS, name = "fc")
fuelA = m.addVars(N, lb = 0.0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = "fuelCaculationA")
m.addConstrs((x[i] + delT * vx[i] == x[i + 1] for i in range(N - 1)), "c1")
m.addConstrs((vx[i] + delT * ux[i] == vx[i + 1] for i in range(N - 1)), "c3")
m.addConstr(x[0] == x0, name="initx")
m.addConstr(vx[0] == vx0, name="initvx")
m.addConstr(ux[0] == ux0, name="initux")
m.addConstr(x[N - 1] == xN, name="finalx")
m.addConstr(vx[N - 1] == vxN, name="finalvx")
for i in range(N -1):
m.addGenConstrPow(vx[i],threeOrder[i],3)
m.addConstr(twoOrder[i] == ux[i] * ux[i])
for i in range(N-1):
m.addConstr(power[i] == d1 * vx[i] + d2 * vx[i] *vx[i] + d3 * threeOrder[i] +
(mass * ux[i]*vx[i])/1000)
m.addConstr(fuelA[i] == alpha + beta1 * power[i] + (beta2 * mass * twoOrder[i] * vx[i]) / 1000)
m.addConstr((fuelBinary[i] == 1 ) >> (fuelConsumption[i] == fuelA[i]))
m.addConstr((fuelBinary[i] == 0) >> (fuelConsumption[i] == alpha))
m.addConstr(power[i] >= 0+0.000001- bigM * (1-fuelBinary[i]))
m.addConstr(power[i] <= 0 + bigM * fuelBinary[i])
expr = LinExpr()
for i in range(N):
expr += fuelConsumption[i]
m.setObjective(expr, GRB.MINIMIZE)
m.optimize()
0
-
-
Gurobi Staff
The following two articles are useful for determining the source of the infeasibility:
- How do I determine why my model is infeasible?
- How do I use 'compute IIS' to find a subset of constraints that are causing model infeasibility?
Following the instructions in the second article, we can compute an Irreducible Inconsistent Subsystem (IIS). This is a subset of constraints and variable bounds that together are infeasible, but become feasible when any single constraint/bound is removed from the subsystem. To do this, call Model.computeIIS() after the call to Model.optimize(), then write out the resulting IIS to disk as an ILP file:
m.optimize()
if m.status == GRB.INFEASIBLE:
m.computeIIS()
m.write("model.ilp")Inspecting the \(\texttt{model.ilp}\) file, we see the following IIS:
\ Model lqrmodel_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
0 power[0]
Subject To
Bounds
power[0] >= inf
EndThis tells us one source of the model's infeasibility: the lower bounds on the \(\texttt{power}\) variables are infinity.
power = m.addVars(N, lb=float('inf'), ub=float("inf"), vtype=GRB.CONTINUOUS, name="power")Even if you fix the lower bounds on the \(\texttt{power}\) variables, the resulting model is still infeasible. After computing another IIS, we see it's not possible to satisfy the \(\texttt{c1}\) and \(\texttt{c3}\) constraints given the initial values \(\texttt{x0}\), \(\texttt{vx0}\), and \(\texttt{ux0}\), the final value \(\texttt{xN}\), and the lower bounds of \(-1\) on the \(\texttt{ux}\) variables:
\ Model lqrmodel_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
Subject To
c1[0]: x[0] - x[1] + vx[0] = 0
c1[1]: x[1] - x[2] + vx[1] = 0
c1[2]: x[2] - x[3] + vx[2] = 0
c1[3]: x[3] - x[4] + vx[3] = 0
c1[4]: x[4] - x[5] + vx[4] = 0
c1[5]: x[5] - x[6] + vx[5] = 0
c1[6]: x[6] - x[7] + vx[6] = 0
c1[7]: x[7] - x[8] + vx[7] = 0
c1[8]: x[8] - x[9] + vx[8] = 0
c1[9]: x[9] - x[10] + vx[9] = 0
c1[10]: x[10] - x[11] + vx[10] = 0
c1[11]: x[11] - x[12] + vx[11] = 0
c1[12]: x[12] - x[13] + vx[12] = 0
c1[13]: x[13] - x[14] + vx[13] = 0
c1[14]: x[14] - x[15] + vx[14] = 0
c1[15]: x[15] - x[16] + vx[15] = 0
c1[16]: x[16] - x[17] + vx[16] = 0
c3[0]: vx[0] - vx[1] + ux[0] = 0
c3[1]: vx[1] - vx[2] + ux[1] = 0
c3[2]: vx[2] - vx[3] + ux[2] = 0
c3[3]: vx[3] - vx[4] + ux[3] = 0
c3[4]: vx[4] - vx[5] + ux[4] = 0
c3[5]: vx[5] - vx[6] + ux[5] = 0
c3[6]: vx[6] - vx[7] + ux[6] = 0
c3[7]: vx[7] - vx[8] + ux[7] = 0
c3[8]: vx[8] - vx[9] + ux[8] = 0
c3[9]: vx[9] - vx[10] + ux[9] = 0
c3[10]: vx[10] - vx[11] + ux[10] = 0
c3[11]: vx[11] - vx[12] + ux[11] = 0
c3[12]: vx[12] - vx[13] + ux[12] = 0
c3[13]: vx[13] - vx[14] + ux[13] = 0
c3[14]: vx[14] - vx[15] + ux[14] = 0
c3[15]: vx[15] - vx[16] + ux[15] = 0
initx: x[0] = 1
initvx: vx[0] = 25
initux: ux[0] = 0.1
finalx: x[17] = 300
Bounds
x[0] free
x[1] free
x[2] free
x[3] free
x[4] free
x[5] free
x[6] free
x[7] free
x[8] free
x[9] free
x[10] free
x[11] free
x[12] free
x[13] free
x[14] free
x[15] free
x[16] free
x[17] free
vx[0] free
vx[1] free
vx[2] free
vx[3] free
vx[4] free
vx[5] free
vx[6] free
vx[7] free
vx[8] free
vx[9] free
vx[10] free
vx[11] free
vx[12] free
vx[13] free
vx[14] free
vx[15] free
vx[16] free
ux[0] free
ux[1] >= -1
ux[2] >= -1
ux[3] >= -1
ux[4] >= -1
ux[5] >= -1
ux[6] >= -1
ux[7] >= -1
ux[8] >= -1
ux[9] >= -1
ux[10] >= -1
ux[11] >= -1
ux[12] >= -1
ux[13] >= -1
ux[14] >= -1
ux[15] >= -1
EndDo the \(\texttt{c1}\) and \(\texttt{c3}\) constraints look correct to you? Are you aware of a feasible solution to this problem? If so, we can identify the constraint(s) violated by the feasible solution. Those constraints would need to be corrected.
0 -
-
Thank you very much Mr. Towle,
I have revised the bounds of power variable. And, according to the code warning, I set the mode.Param.NonCovex=2. Now, the model does not give me the model infeasible error. But the model is running forever without giving the answer. Here is the updated code below. Thanks.
def OptimalTrajectory1(timesteps=20, initialSpeed=15, finalSpeed = 15, distance=300, minimumSpeed=1, maximumSpeed=30):
### Zhihong Yao integrated ....
# Parameters (constants)
N = timesteps # Instants or Steps
delT = 1 # Time duration
u_max_x = 2
# state and input variables initialization (variables)
x0 = 0
vx0 = initialSpeed
ux0 = 0
uxN = 0
# final state to reach at:
xN = distance
vxN = finalSpeed
solutions = {}
alpha = 0.666
beta1 = 0.072
beta2 = 0.0344
d1 = 0.269
d2 = 0.0171
d3 = 0.000672
mass = 1680
bigM = 5000
# Create a new model
m = Model("lqrmodel")
m.Params.OutputFlag = 0
m.Params.NonConvex = 2
m.Params.TimeLimit =10
# Create variables
x = m.addVars(N, lb=0.0, ub=float('inf'), vtype=GRB.CONTINUOUS, name="x")
vx = m.addVars(N, lb=minimumSpeed, ub=maximumSpeed, vtype=GRB.CONTINUOUS, name="vx")
ux = m.addVars(N, lb=-u_max_x, ub=u_max_x, vtype=GRB.CONTINUOUS, name="ux")
fuelBinary = m.addVars(N, vtype = GRB.BINARY, name="FCIndicator")
power = m.addVars(N, lb=0, ub=float("inf"), vtype=GRB.CONTINUOUS, name="power")
threeOrder = m.addVars(N, lb = 0, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
twoOrder = m.addVars(N, lb = 0, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
fuelConsumption = m.addVars(N, lb= 0.00, ub = float('inf'), vtype = GRB.CONTINUOUS, name = "fc")
fuelA = m.addVars(N, lb = 0.0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = "fuelCaculationA")
m.addConstrs((x[i] + delT * vx[i] == x[i + 1] for i in range(N - 1)), "c1")
m.addConstrs((vx[i] + delT * ux[i] == vx[i + 1] for i in range(N - 1)), "c3")
m.addConstr(x[0] == x0, name="initx")
m.addConstr(vx[0] == vx0, name="initvx")
m.addConstr(ux[0] == ux0, name="initux")
m.addConstr(ux[N-1]== uxN, name = "finalux")
m.addConstr(x [N - 1] == xN, name="finalx")
m.addConstr(vx[N - 1] == vxN, name="finalvx")
for i in range(N -1):
m.addGenConstrPow(vx[i],threeOrder[i],3)
m.addConstr(twoOrder[i] == ux[i] * ux[i])
for i in range(N-1):
m.addConstr(power[i] == d1 * vx[i] + d2 * vx[i] *vx[i] + d3 * threeOrder[i] +
(mass * ux[i]*vx[i])/1000)
m.addConstr(fuelA[i] == alpha + beta1 * power[i] + (beta2 * mass * twoOrder[i] * vx[i]) / 1000)
m.addConstr((fuelBinary[i] == 1 ) >> (fuelConsumption[i] == fuelA[i]))
m.addConstr((fuelBinary[i] == 0) >> (fuelConsumption[i] == alpha))
m.addConstr(power[i] >= 0+0.000001- bigM * (1-fuelBinary[i]))
m.addConstr(power[i] <= 0 + bigM * fuelBinary[i])
expr = LinExpr()
for i in range(N):
expr += fuelConsumption[i]
m.setObjective(expr, GRB.MINIMIZE)
m.optimize()
if m.status == GRB.INFEASIBLE:
m.computeIIS()
m.write("model.ilp")
print(m.status)
if m.status == 2:
# figure = plt.figure()
# ax = figure.add_subplot(1,2,1)
# ax1 = figure.add_subplot(1,2,2)
plt.plot([v.X for i, v in vx.items()])
solutions[vxN] = [v.x for i, v in vx.items()]
print(vxN, "has optimal results")
plt.show()
# plt.plot([xn.X for i,xn in x.items()])
return solutions0 -
-
-
You set the TimeLimit parameter to 10, so Gurobi should not spend more than 10 seconds solving the problem. Please post Gurobi's output, which is displayed if you leave the OutputFlag parameter at its default value. Also, what arguments are you passing to the \(\texttt{OptimalTrajectory1}\) the function?
0 -
-
This is the output of the model. It takes hours to obtain a proper solution. The arguments I passed to the OptimalTrajectory1 is the default values defined in the definition of the function, i.e.,
timesteps=20, initialSpeed=15, finalSpeed = 15, distance=300, minimumSpeed=1, maximumSpeed=30. The time to find a solution for this model is unexpected large, I suspect that I must did something wrong~
Optimize a model with 82 rows, 180 columns and 196 nonzeros
Model fingerprint: 0xead3ea56
Model has 57 quadratic constraints
Model has 57 general constraints
Variable types: 160 continuous, 20 integer (20 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+03]
QMatrix range [2e-02, 2e+00]
QLMatrix range [7e-04, 1e+00]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 3e+01]
RHS range [2e+01, 5e+03]
QRHS range [7e-01, 7e-01]
GenCon rhs range [7e-01, 7e-01]
GenCon coe range [1e+00, 1e+00]
Presolve added 132 rows and 2958 columns
Presolve time: 0.02s
Presolved: 472 rows, 3172 columns, 9888 nonzeros
Presolved model has 17 SOS constraint(s)
Presolved model has 69 bilinear constraint(s)
Variable types: 3100 continuous, 72 integer (72 binary)
Root relaxation: objective 1.338484e+01, 239 iterations, 0.01 seconds (0.01 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 13.38484 0 75 - 13.38484 - - 0s
0 0 13.38484 0 81 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 84 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 103 - 13.38484 - - 0s
0 0 13.38484 0 102 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 98 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 97 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 89 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 104 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 101 - 13.38484 - - 0s
0 0 13.38484 0 101 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 97 - 13.38484 - - 0s
0 0 14.08274 0 130 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 16.98143 0 145 - 16.98143 - - 0s
0 0 16.98143 0 147 - 16.98143 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.96743 0 146 - 18.96743 - - 0s
0 0 19.30901 0 143 - 19.30901 - - 0s
0 0 19.31416 0 154 - 19.31416 - - 0s
0 0 19.58767 0 149 - 19.58767 - - 0s
0 0 19.59163 0 150 - 19.59163 - - 0s
0 0 19.93134 0 149 - 19.93134 - - 0s
0 0 19.94079 0 149 - 19.94079 - - 0s
0 0 19.94594 0 152 - 19.94594 - - 0s
0 0 19.97007 0 151 - 19.97007 - - 0s
H 0 0 28.9742352 19.97007 31.1% - 0s
0 0 19.97007 0 84 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s0 -
-
-
Is there more to the log file? This log snippet shows less than one second's worth of Gurobi's solve. The log is also missing some important information at the beginning, like Gurobi version, operating system, and which parameters are set.
Gurobi finds an initial feasible solution with objective value 28.9742352 very quickly, but I can't tell what happens afterwards without the rest of the log file. At any rate, Gurobi should not spend more than 10 seconds solving, because you set the TimeLimit parameter to 10.
0 -
-
Hi Towle,
I set the MIPGap parameter to see whether I can obtain a proper solution, instead of the time limit. Here is the updated code.
def OptimalTrajectory1(timesteps=20, initialSpeed=15, finalSpeed = 15, distance=300, minimumSpeed=1, maximumSpeed=30):
### Zhihong Yao integrated ....
# Parameters (constants)
N = timesteps # Instants or Steps
delT = 1 # Time duration
u_max_x = 2
# state and input variables initialization (variables)
x0 = 0
vx0 = initialSpeed
ux0 = 0
uxN = 0
# final state to reach at:
xN = distance
vxN = finalSpeed
solutions = {}
alpha = 0.666
beta1 = 0.072
beta2 = 0.0344
d1 = 0.269
d2 = 0.0171
d3 = 0.000672
mass = 1680
bigM = 5000
# Create a new model
m = Model("lqrmodel")
# m.Params.OutputFlag = 0
m.Params.NonConvex = 2
m.Params.MIPGap =0.03
# m.Params.TimeLimit =10
# Create variables
x = m.addVars(N, lb=0.0, ub=float('inf'), vtype=GRB.CONTINUOUS, name="x")
vx = m.addVars(N, lb=minimumSpeed, ub=maximumSpeed, vtype=GRB.CONTINUOUS, name="vx")
ux = m.addVars(N, lb=-u_max_x, ub=u_max_x, vtype=GRB.CONTINUOUS, name="ux")
fuelBinary = m.addVars(N, vtype = GRB.BINARY, name="FCIndicator")
power = m.addVars(N, lb=0, ub=float("inf"), vtype=GRB.CONTINUOUS, name="power")
threeOrder = m.addVars(N, lb = 0, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
twoOrder = m.addVars(N, lb = 0, ub = GRB.INFINITY,vtype=GRB.CONTINUOUS)
fuelConsumption = m.addVars(N, lb= 0.00, ub = float('inf'), vtype = GRB.CONTINUOUS, name = "fc")
fuelA = m.addVars(N, lb = 0.0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = "fuelCaculationA")
m.addConstrs((x[i] + delT * vx[i] == x[i + 1] for i in range(N - 1)), "c1")
m.addConstrs((vx[i] + delT * ux[i] == vx[i + 1] for i in range(N - 1)), "c3")
m.addConstr(x[0] == x0, name="initx")
m.addConstr(vx[0] == vx0, name="initvx")
m.addConstr(ux[0] == ux0, name="initux")
m.addConstr(ux[N-1]== uxN, name = "finalux")
m.addConstr(x [N - 1] == xN, name="finalx")
m.addConstr(vx[N - 1] == vxN, name="finalvx")
for i in range(N -1):
m.addGenConstrPow(vx[i],threeOrder[i],3)
m.addConstr(twoOrder[i] == ux[i] * ux[i])
for i in range(N-1):
m.addConstr(power[i] == d1 * vx[i] + d2 * vx[i] *vx[i] + d3 * threeOrder[i] +
(mass * ux[i]*vx[i])/1000)
m.addConstr(fuelA[i] == alpha + beta1 * power[i] + (beta2 * mass * twoOrder[i] * vx[i]) / 1000)
m.addConstr((fuelBinary[i] == 1 ) >> (fuelConsumption[i] == fuelA[i]))
m.addConstr((fuelBinary[i] == 0) >> (fuelConsumption[i] == alpha))
m.addConstr(power[i] >= 0+0.000001- bigM * (1-fuelBinary[i]))
m.addConstr(power[i] <= 0 + bigM * fuelBinary[i])
expr = LinExpr()
for i in range(N):
expr += fuelConsumption[i]
m.setObjective(expr, GRB.MINIMIZE)
m.optimize()
if m.status == GRB.INFEASIBLE:
m.computeIIS()
m.write("model.ilp")
print(m.status)
if m.status == 2:
# figure = plt.figure()
# ax = figure.add_subplot(1,2,1)
# ax1 = figure.add_subplot(1,2,2)
plt.plot([v.X for i, v in vx.items()])
solutions[vxN] = [v.x for i, v in vx.items()]
print(vxN, "has optimal results")
plt.show()
# plt.plot([xn.X for i,xn in x.items()])
return solutionsand here is the full output
Restricted license - for non-production use only - expires 2024-10-28
Set parameter NonConvex to value 2
Set parameter MIPGap to value 0.03
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64)
CPU model: AMD Ryzen 5 4600H with Radeon Graphics, instruction set [SSE2|AVX|AVX2]
Thread count: 6 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 82 rows, 180 columns and 196 nonzeros
Model fingerprint: 0xead3ea56
Model has 57 quadratic constraints
Model has 57 general constraints
Variable types: 160 continuous, 20 integer (20 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+03]
QMatrix range [2e-02, 2e+00]
QLMatrix range [7e-04, 1e+00]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 3e+01]
RHS range [2e+01, 5e+03]
QRHS range [7e-01, 7e-01]
GenCon rhs range [7e-01, 7e-01]
GenCon coe range [1e+00, 1e+00]
Presolve added 132 rows and 2958 columns
Presolve time: 0.01s
Presolved: 472 rows, 3172 columns, 9888 nonzeros
Presolved model has 17 SOS constraint(s)
Presolved model has 69 bilinear constraint(s)
Variable types: 3100 continuous, 72 integer (72 binary)
Root relaxation: objective 1.338484e+01, 239 iterations, 0.01 seconds (0.01 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 13.38484 0 75 - 13.38484 - - 0s
0 0 13.38484 0 81 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 84 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 79 - 13.38484 - - 0s
0 0 13.38484 0 103 - 13.38484 - - 0s
0 0 13.38484 0 102 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 98 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 97 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 96 - 13.38484 - - 0s
0 0 13.38484 0 95 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 90 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 89 - 13.38484 - - 0s
0 0 13.38484 0 91 - 13.38484 - - 0s
0 0 13.38484 0 104 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 101 - 13.38484 - - 0s
0 0 13.38484 0 101 - 13.38484 - - 0s
0 0 13.38484 0 107 - 13.38484 - - 0s
0 0 13.38484 0 97 - 13.38484 - - 0s
0 0 14.08274 0 130 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 14.08274 0 131 - 14.08274 - - 0s
0 0 16.98143 0 145 - 16.98143 - - 0s
0 0 16.98143 0 147 - 16.98143 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 17.13380 0 144 - 17.13380 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.84705 0 146 - 18.84705 - - 0s
0 0 18.96743 0 146 - 18.96743 - - 0s
0 0 19.30901 0 143 - 19.30901 - - 0s
0 0 19.31416 0 154 - 19.31416 - - 0s
0 0 19.58767 0 149 - 19.58767 - - 0s
0 0 19.59163 0 150 - 19.59163 - - 0s
0 0 19.93134 0 149 - 19.93134 - - 0s
0 0 19.94079 0 149 - 19.94079 - - 0s
0 0 19.94594 0 152 - 19.94594 - - 0s
0 0 19.97007 0 151 - 19.97007 - - 0s
H 0 0 28.9742352 19.97007 31.1% - 0s
0 0 19.97007 0 84 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 19.97007 0 81 28.97424 19.97007 31.1% - 0s
0 0 21.72067 0 83 28.97424 21.72067 25.0% - 0s
0 0 23.51405 0 89 28.97424 23.51405 18.8% - 1s
0 0 23.85033 0 87 28.97424 23.85033 17.7% - 1s
0 0 23.89591 0 88 28.97424 23.89591 17.5% - 1s
0 0 23.95859 0 87 28.97424 23.95859 17.3% - 1s
0 0 23.96066 0 94 28.97424 23.96066 17.3% - 1s
0 0 24.02787 0 96 28.97424 24.02787 17.1% - 1s
0 0 24.03217 0 96 28.97424 24.03217 17.1% - 1s
0 0 24.06199 0 95 28.97424 24.06199 17.0% - 1s
0 0 24.06204 0 95 28.97424 24.06204 17.0% - 1s
0 0 24.06453 0 94 28.97424 24.06453 16.9% - 1s
0 0 24.06775 0 94 28.97424 24.06775 16.9% - 1s
0 0 24.06775 0 96 28.97424 24.06775 16.9% - 1s
0 0 24.09530 0 85 28.97424 24.09530 16.8% - 1s
H 0 0 28.6004176 24.09530 15.8% - 1s
0 2 24.09530 0 85 28.60042 24.09530 15.8% - 1s
* 1128 876 147 28.4301447 24.11471 15.2% 24.1 2s
H 1318 894 28.2074716 25.84641 8.37% 26.4 5s
1400 952 26.84018 33 91 28.20747 26.84018 4.85% 32.6 10s
* 2927 1567 208 28.1034988 27.06250 3.70% 35.1 13s
4234 2075 27.27980 62 126 28.10350 27.21741 3.15% 33.4 15s
Cutting planes:
Gomory: 19
Implied bound: 14
Projected implied bound: 15
MIR: 35
Flow cover: 13
RLT: 59
Relax-and-lift: 37
PSD: 6
Explored 6213 nodes (205116 simplex iterations) in 16.66 seconds (16.63 work units)
Thread count was 12 (of 12 available processors)
Solution count 5: 28.1035 28.2075 28.4301 ... 28.9742
Optimal solution found (tolerance 3.00e-02)
Warning: max constraint violation (1.8618e-04) exceeds tolerance
Warning: max general constraint violation (1.8618e-04) exceeds tolerance
Best objective 2.810349875926e+01, best bound 2.726255231969e+01, gap 2.9923%
2
15 has optimal results
Process finished with exit code 00 -
-
-
Gurobi solves the model to a 3% optimality gap in 16.66 seconds. Are you having any other issues?
I see you're using Model.addGenConstrPow() to add the constraints \(\texttt{threeOrder[i]} = \texttt{vx[i]}^3\). With this approach, Gurobi builds a piecewise-linear approximation of the function \(x^3\). Alternatively, you could try modeling this cubic relationship directly by introducing two additional non-convex constraints. In particular, comment out the Model.addGenConstrPow() line. Add a set of auxiliary non-negative variables \(\texttt{aux}\):
aux = m.addVars(N, ub=maximumSpeed**2)
Then, model \(\texttt{threeOrder[i]} = \texttt{vx[i]}^3\) using two bilinear constraints and these new auxiliary variables:
for i in range(N -1):
m.addConstr(aux[i] == vx[i] * vx[i])
m.addConstr(threeOrder[i] == aux[i] * vx[i])You could check if the solutions/performance are better with this approach.
It's good practice to provide the tightest possible bounds for variables participating in bilinear constraints like this. See the article How do I model multilinear terms in Gurobi? for more details.
0 -
-
Thank you for your professional input. I will continue to work on this problem. I currently have one more issue, for which I made a new post in this platform. Thanks again for your kind help.
0 -
Please sign in to leave a comment.
Comments
8 comments