Need guidance in finding error in transport optimization model

I'm writing a H2 transportation & production cost optimization model and it keeps returning the following:

Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64) 
 
CPU model: 12th Gen Intel(R) Core(TM) i9-12900K, instruction set [SSE2|AVX|AVX2] 
Thread count: 16 physical cores, 24 logical processors, using up to 24 threads  

Optimize a model with 96006 rows, 42492 columns and 1336724 nonzeros 
Coefficient statistics:
  Matrix range [1e+00, 1e+00]
  Objective range [2e+06, 1e+10]
  Bounds range [0e+00, 0e+00]
  RHS range [9e-02, 1e+10] 
Warning: Model contains large objective coefficients 
Warning: Model contains large rhs
  Consider reformulating model or setting NumericFocus parameter to avoid numerical issues.
 Iteration   Objective       Primal Inf.   Dual Inf. Time
  0         4.5300787e+15   4.403554e+07   0.000000e+00 0s

  Solved in 0 iterations and 0.05 seconds (0.03 work units)
 Infeasible model

#Gurobi Variables---------------------------------------------------
model = Model("Japan_2030")

Product_PV = {} #Production of Solar PV power derived H2
Product_ON = {} #Production of Onshore wind power derived H2
Product_OFF = {} #Production of Offshore wind power derived H2
Product_IM = {} #Quantity of Imported H2
Traffic_PV = {} #Traffic of Solar PV power derived H2
Traffic_ON = {} #Traffic of Onshore wind power derived H2
Traffic_OFF = {} #Traffic of Offshore wind power derived H2
Traffic_IM = {} #Traffic of Imported H2

for i in tqdm(I):
    Product_PV[i] = model.addVar(vtype="C", name="Product_PV(%s)"%(i))
    Product_ON[i] = model.addVar(vtype="C", name="Product_ON(%s)"%(i))
    Product_OFF[i] = model.addVar(vtype="C", name="Product_OFF(%s)"%(i))
    Product_IM[i] = model.addVar(vtype="C", name="Product_IM(%s)"%(i))

    for j in J:
        if AdjData.loc[i,j] == 1:
            Traffic_PV[i,j] = model.addVar(vtype="C", name="Traffic_PV(%s,%s)"%(i,j))
            Traffic_ON[i,j] = model.addVar(vtype="C", name="Traffic_ON(%s,%s)"%(i,j))
            Traffic_OFF[i,j] = model.addVar(vtype="C", name="Traffic_OFF(%s,%s)"%(i,j))
            Traffic_IM[i,j] = model.addVar(vtype="C", name="Traffic_IM(%s,%s)"%(i,j))

model.update()


#Gurobi Constraints---------------------------------------------------
for i in tqdm(I):
    model.addConstr(Product_PV[i] + Product_ON[i] + Product_OFF[i] + Product_IM[i]
                      + quicksum(Traffic_PV[i,j] + Traffic_ON[i,j] + Traffic_OFF[i,j] + Traffic_IM[i,j]
                                  - Traffic_PV[j,i] - Traffic_ON[j,i] - Traffic_OFF[j,i] - Traffic_IM[j,i]
                                   if AdjData.loc[i,j] == 1 else 0 for j in J) == DEMAND[i], name="Demand(%s)"%i)

for j in tqdm(J):
    model.addConstr(Product_PV[j] + quicksum(Traffic_PV[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_PV[j], name="Capacity_PV(%s)"%j)
    model.addConstr(Product_ON[j] + quicksum(Traffic_ON[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_ON[j], name="Capacity_ON(%s)"%j)
    model.addConstr(Product_OFF[j] + quicksum(Traffic_OFF[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_OFF[j], name="Capacity_OFF(%s)"%j)
    model.addConstr(Product_IM[j] + quicksum(Traffic_IM[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_IM[j], name="Capacity_IM(%s)"%j)

#Gurobi Objective---------------------------------------------------
model.setObjective(quicksum(LCOH_PV_LIST[j] * Product_PV[j] + LCOH_ON_LIST[j] * Product_ON[j]
                            + LCOH_OFF_LIST[j] * Product_OFF[j] + LCOH_IM_LIST[j] * Product_IM[j] for j in J)
                    +quicksum( TRANSPORTCOSTLIST[i,j] * (Traffic_PV[i,j] + Traffic_ON[i,j] + Traffic_OFF[i,j] + Traffic_IM[i,j])
                                if AdjData.loc[i,j] == 1 else 0 for (i,j) in tuple_list), GRB.MINIMIZE)

model.optimize()

I've limited Traffic_ variables to only adjacent municipals(I and J are both city codes, tuple_list is the tuple version). AdjData is municipal adjacency data (1 if adjacent, else 0). LCOH_ lists and TRANSPORTCOSTLIST are costs for producing and transporting H2 (both dict). 

Another algorithm where I didn't consider adjacency of municipals gave me a realistic optimal value. Is there something I'm doing wrong?

UPDATE: I found I can bypass the error when I omit any one of these 4 contraints. Can having these 4 constraints somehow cause any trouble?

for j in tqdm(J):
  model.addConstr(Product_PV[j] + quicksum(Traffic_PV[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_PV[j], name="Capacity_PV(%s)"%j)
  model.addConstr(Product_ON[j] + quicksum(Traffic_ON[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_ON[j], name="Capacity_ON(%s)"%j)
  model.addConstr(Product_OFF[j] + quicksum(Traffic_OFF[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_OFF[j], name="Capacity_OFF(%s)"%j)
  model.addConstr(Product_IM[j] + quicksum(Traffic_IM[i,j] if AdjData.loc[i,j] == 1 else 0 for i in I)  <= M_IM[j], name="Capacity_IM(%s)"%j)