Gurobi says model is infeasible, yet I am able to find a solution by hand

Question

I have the following model, which Gurobi says is infeasible, yet I am able to find a solution by hand. Please note that for num_steps of 49, the problem does have a solution.

from gurobipy import Model, GRB

f = Model()

num_steps = 50
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(-float('inf'), 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)
    
K = {}
for step in steps_for_C:
    K[step] = f.addVar(0, float('inf'), name="K[%i]" % step)

N = {}
for step in steps_for_C:
    N[step] = f.addVar(0, float('inf'), name="N[%i]" % step)
    
M = {}
for step in steps_for_C:
    M[step] = f.addVar(0, float('inf'), name="M[%i]" % step)
    
for step in steps:
    L = 10
    f.addConstr(D[0] == 1, name = 'D == 1 initial')
    f.addConstr(C[0] == 0.8*A, name = 'C initial')
    f.addConstr(K[0] == 0, name = 'K initial')
    
    f.addConstr(C[step + 1] == C[step] * 0.9999916666666666 + 0.95*E[step] - F[step]/0.95, 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(N[step] == 0.75 * J[step] - 0.02 * C[step] + 0.02, name = 'N constr')
    f.addGenConstrAbs(M[step], N[step], name = 'M constr')
    f.addConstr(K[step+1] == K[step] + (0.5*M[step]/4500), name = 'K constr')
    f.addConstr(D[step + 1] == D[step] - (K[step] + (step)/113880)*0.2, name = 'D constr')
    f.addConstr(C[step] <= 0.8*A*D[step], 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] <= 0.8*A*D[step] - C[step], name = 'E limit')
    f.addConstr(F[step] <= B, name = 'F less than B')
    f.addConstr(F[step] <= 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()

print()
print('A', A.x)
print('B', B.x)

Here's the solution that I seem to be able to find by hand:

A = 500
B = 2

D = 1
C = 0.8*A

F = B
E = 0

K = 0

for step in steps:
    J = abs(E - F)
    C = C * 0.9999916666666666 + 0.95*E - F/0.95
    K = K + 0.5*abs(0.75 * J - 0.02 * C + 0.02)/4500
    D = D - (K + (step+1)/113880)*0.2
    G = E - F + L
    
    print('G is: ', G)
    print('D is: ', D)
    print('C is: ', C)
    print('----------')

Am I missing something? Why is the model infeasible for 50 timesteps but is feasible for 49? Any help is appreciated.