Gurobi says the model is infeasible, yet I can come up with a solution by hand. Python implementation

Question

Here is my model. Gurobi says its infeasible. My apologies if its hard to follow.

from gurobipy import Model, GRB

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, 8, name="G[%i]" % step)

for step in steps:
    L = 10
    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 + step/1000), 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

L = 10

for step in [1,2,3]:
    J = abs(E - F)
    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. 

Answer 1

Turning the bounds of H from f.addVar(0, float('inf'), name="H[%i]" % step) to f.addVar(-float('inf'), float('inf'), name="H[%i]" % step) fixes this.