Gurobi says the model is infeasible, yet I can come up with a solution by hand. Python implementation
Answered
Here is my model. Gurobi says its infeasible. My apologies if its hard to follow.
f = Model()
num_steps = 3
steps = [*range(num_steps)]
steps_for_C = [*range(num_steps + 1)]
A = f.addVar(0, float('inf'), name = 'A')
B = f.addVar(0, float('inf'), name = 'B')
C = {}
for step in steps_for_C:
C[step] = f.addVar(0, float('inf'), name="C[%i]" % step)
D = {}
for step in steps_for_C:
D[step] = f.addVar(0, 1, name="D[%i]" % step)
E = {}
for step in steps_for_C:
E[step] = f.addVar(0, float('inf'), name="E[%i]" % step)
F = {}
for step in steps_for_C:
F[step] = f.addVar(0, float('inf'), name="F[%i]" % step)
H = {}
for step in steps_for_C:
H[step] = f.addVar(0, float('inf'), name="H[%i]" % step)
J = {}
for step in steps_for_C:
J[step] = f.addVar(0, float('inf'), name="J[%i]" % step)
G = {}
for step in steps:
G[step] = f.addVar(0, 80, name="G[%i]" % step)
for step in steps:
L = 100
f.addConstr(D[0] == 1., name = 'C == 1')
f.addConstr(C[0] == 0.8*A, name = 'C initial')
f.addConstr(C[step + 1] == C[step] * 0.9999 + E[step+1] - F[step+1], name = 'C update')
f.addConstr(H[step] == E[step] - F[step], name = 'H constr')
f.addGenConstrAbs(J[step], H[step], name = 'J constr')
f.addConstr(D[step + 1] == D[step] - (0.75 * J[step+1] - 0.02 * C[step+1] + 0.02) * 0.001, name = 'D constr')
f.addConstr(C[step] <= 0.8*A*D[step+1], name = 'C max')
f.addConstr(C[step] >= 0.2*A, name = 'C min')
f.addConstr(E[step] <= B, name = 'E less than B')
f.addConstr(E[step+1] <= 0.8*A*D[step+1] - C[step], name = 'E limit')
f.addConstr(F[step] <= B, name = 'F less than B')
f.addConstr(F[step+1] <= C[step] - 0.2*A, name = 'F limit')
f.addConstr(-L + G[step] - E[step] + F[step] == 0, name = 'p constraint')
f.setObjective(A, GRB.MINIMIZE)
f.params.FeasibilityTol = 0.01
f.params.NonConvex = 2
f.optimize()However, I can easily come up with a solution by hand, as in:
A = 10
B = 2
D = 1
C = 8
F = 2
E = 0
for step in [1,2,3]:
J = abs(PB_charging - PB_discharging)
D -= ((0.75 * J - 0.02 * C + 0.02) * 0.001 + step/1000)
C = C * 0.9999 + E - F
G = E - F + L
print('G is: ', G)
print('D is: ', D)
print('C is: ', C)
print('----------')Since I can come up with a solution by hand, why is the model infeasible? Am I missing something here? Perhaps the order in which I arrange my constraints is wrong? Any help is appreciated.
-
For easy-to-follow code, please see here: https://stackoverflow.com/questions/69600967/gurobi-says-the-model-is-infeasible-yet-i-can-come-up-with-a-solution-by-hand
0 -
Hi Vadim,
The solution you compute by hand has
G[1] is: 98
However, you set the upper bound of variable \(\texttt{G}\) to \(80\) so the solution you computed by hand is infeasible. Indeed, setting the upper bound of \(\texttt{G}\) to \(100\) makes the problem feasible.
In your other post you stated that you cannot compute the IIS. You have to compute the IIS after your problem has been declared to be infeasible, i.e.,
f.optimize()
f.computeIIS()The computed IIS is then
Subject To
H_constr: - E[2] + F[2] + H[2] = 0
p_constraint: - E[2] + F[2] + G[2] = 100
Bounds
E[2] free
F[2] free
-infinity <= G[2] <= 80
Endwhich points to some issue with variable bounds.
Best regards,
Jaromił0 -
Thank you very much Jaromil. My mistake, I should've paid more attention to my code.
However, if I change my L to 10, and my G upper limit to 8, I still do get a solution by hand. Is there anything that you believe I've done wrong?
0 -
When I make these changes I still do get the 'model infeasible' error.
0 -
When I run the code to determine constrains that aren't satisfied, I only get the 'H constraint' that does not work:
print('\nThe following constraint cannot be satisfied:')
removed = []
for c in f.getConstrs():if c.IISConstr:
print('%s' % c.constrName)But how can 'H constraint' not be satisfied? Its simple the absolute difference betwween the E and F.
0 -
I should mention, in my simulation, some of the model variables (like J) are calculated based on other variables that I am actually trying to optimize. J being the absolute difference between E and F should be simply calculated, not optimized. I do not know of any other way to do it in Gurobi, other than to make it into a variable. Please let me know if there is a simpler way to managing something like this.
0 -
Hi Vadim,
When you set \(G=8, E=0, F=2\), you then get
\[\begin{align}
- 0 + 2 + H &= 0\\
- 0 + 2 + 8 &= 10
\end{align}\]so \(H=-2\). But you define \(H\geq 0\). Setting the lower bound of \(H\) to \(- \infty\) should solve the problem.
Best regards,
Jaromił0 -
Thank you Jaromil. That solved the previous issue.
Now I do have another question. When I increase the number of steps (num_steps) from 3 to 12 I get a solution. When I increase it past 12 to 13 Gurobi says the model is infeasible and it appears multiple constraints are violated.
Is there a way to see by how much each of those constraints gets violated? In my setup I do not believe the problem is unsolvable, I think the issue is numerical (precision-related). I need my answers accurate only to 2 or 3 decimal points, except for variable D. D is usually "smooth" and needs to continually decrease and never increase.
0 -
Additionally, I should be able to pick an A large enough so that it will never become less than 0 and it appears this problem should be feasible.
0 -
Hi Vadim,
I am not sure if I use your latest version of the code. I am able to set \(\texttt{num_steps}\) up to 58 and still get a feasible solution. \(\texttt{num_steps=59}\) results in an infeasible model. Working with so many steps is very tedious. Could you provide the latest version of your code which turns infeasible if \(\texttt{num_steps=13}\). Did you have a look at the feasRelax function? It is a different way of dealing with infeasible models. You can find more info in the Knowledge Base article How do I determine why my model is infeasible?
There is a constraint called \(\texttt{C==1}\) but it states
D[0] == 1
Is this on purpose?
Best regards,
Jaromił0 -
Thank you Jaromil. Yes, absolutely, the latest version is here: https://stackoverflow.com/questions/69640879/gurobi-says-model-is-infeasible-yet-i-am-able-to-find-a-solution-by-hand . It works for 49 steps but not 50. I have fixed the naming for D[0] == 1 constraint. I have tried using feasRelaxS with limited success. On my main model it either takes too long to run or gives me weird answers, like A = 0 and B > 0.
As a side question, if I do use feasRelaxS, how do I access the values of the variables I am trying to optimize? I noticed that when I run f.feasRelaxS(...) followed by f.optimize(), I still can't access variable values like A.x and B.x. Thank you for your help.
0 -
Hi Vadim,
As the user sascha on stackoverflow proposed, you should try to fix the variable bounds to the feasible point you computed and then check why it is infeasible via either IIS computation or the feasRelax method. The fixing can be done via
A.lb = 500
A.ub = 500
B.lb = 2
B.ub = 2
f.update()
D[0].lb = 1
D[0].ub = 1
C[0].lb = 0.8*A.lb
C[0].ub = 0.8*A.lb
K[0].lb = 0
K[0].ub = 0
for step in steps:
f.update()
F[step].lb = B.lb
F[step].ub = B.lb
E[step].lb = 0
E[step].ub = 0
J[step].lb = abs(E[step].lb - F[step].lb)
C[step+1].lb = C[step].lb * 0.9999916666666666 + 0.95*E[step].lb - F[step].lb/0.95
K[step+1].lb = K[step].lb + 0.5*abs(0.75 * J[step].lb - 0.02 * C[step].lb + 0.02)/4500
D[step+1].lb = D[step].lb - (K[step].lb + (step+1)/113880)*0.2
G[step].lb = E[step].lb - F[step].lb + L
J[step].ub = abs(E[step].lb - F[step].lb)
C[step+1].ub = C[step].lb * 0.9999916666666666 + 0.95*E[step].lb - F[step].lb/0.95
K[step+1].ub = K[step].lb + 0.5*abs(0.75 * J[step].lb - 0.02 * C[step].lb + 0.02)/4500
D[step+1].ub = D[step].lb - (K[step].lb + (step+1)/113880)*0.2
G[step].ub = E[step].lb - F[step].lb + L
f.params.FeasibilityTol = 1e-2
f.params.NonConvex = 2
f.optimize()
if f.Status == GRB.INFEASIBLE:
f.computeIIS()
f.write("myIIS.ilp")
f.feasRelaxS(1,False,True,True)
f.optimize()
print('A', A.x)
print('B', B.x)You should definitely double check my code before just using it.
The IIS I got is
Minimize
Subject To
C_update: - 0.9999916666666666 C[0] + C[1] - 0.95 E[0]
+ 1.052631578947368 F[0] = 0
Bounds
-infinity <= C[0] <= 400
C[1] >= 399.9966666666667
-infinity <= E[0] <= 0
F[0] >= 2
Endwhich states that
C[1] = 0.9999916666666666 C[0] - 1.052631578947368 F[0]
is infeasible for your fixed feasible point because
0.9999916666666666 * 400 - 1.052631578947368 * 2 = 397.8914 < 399.996
As a side question, if I do use feasRelaxS, how do I access the values of the variables I am trying to optimize? I noticed that when I run f.feasRelaxS(...) followed by f.optimize(), I still can't access variable values like A.x and B.x
I was not able to reproduce this issue. After calling
f.feasRelaxS(1,False,True,True)
f.optimize()I was able to execute
print('A', A.x)
print('B', B.x)without errors.
Best regards,
Jaromił0 -
Hello Jaromil. Thank you very much for the suggestions. The problem ended up being consistently subtracting (step+1)/113880 term from D which makes it less than 0.
I had another question. Is there a way to tell if the optimization is progressing and the expected end time of the optimization? I'm looking for a way to gather this information after the optimization started.
As an extension to this question, how would I be able to make and to view the gurobi.log file while working in a jupyter notebook? What sort of information can gurobi.log file tell me to understand the progress of my optimization?
0 -
Hi Vadim,
Is there a way to tell if the optimization is progressing and the expected end time of the optimization? I'm looking for a way to gather this information after the optimization started.
In general, there is no way of telling how long an optimization run will take. If you have strict time restrictions, you should always set the TimeLimit parameter. During the optimization, you can track the value of the MIPGap in the output to see how close the solver got to the desired MIPGap. However, this value can be misleading because even if the reported MIPGap value is close to the final desired termination value, it does not mean that the optimization is going to finish soon. It is often the most difficult part to close the final bit of the MIPGap and terminate.
As an extension to this question, how would I be able to make and to view the gurobi.log file while working in a jupyter notebook? What sort of information can gurobi.log file tell me to understand the progress of my optimization?
You can just read the output in the console. If you want to perform automatic evaluation based on the current log output, you can use callbacks to get certain information about the optimization process. You could either read each log line in the MESSAGE callback or get the best bound and best feasible solution value in the MIP callback. You then have to compute the MIPGap on your own within the callback.
Best regards,
Jaromił0
Please sign in to leave a comment.
Comments
14 comments