Too much time wasted in adding quadratic constraint for a quadratic minimization problem

"Gurobipy model formulation"

SEVRP = Model()

SEVRP.setParam('OutputFlag', 0)
SEVRP.setParam('TimeLimit', tlimit)
SEVRP.setParam('NumericFocus', 1)

# 'Auxiliary variables'
# Z = SEVRP.addMVar((num_nodes, num_nodes), vtype=GRB.BINARY, name="Z")
Z = SEVRP.addVars(Z_set, vtype=GRB.BINARY, name="Z")

# Dual vars of DRO
gamma = SEVRP.addVars(D, name="gamma")
alpha = SEVRP.addVars(Z_set, D, name="alpha")
Q = SEVRP.addVars(Z_set, D, lb=0, name="Q")
chi = SEVRP.addVars(k_ad, lb=0, name="chi")
x_kad = SEVRP.addVars(num_x[0], D, k_ad, vtype=GRB.BINARY, name="x_kad")
zeta_s = SEVRP.addVars(A.shape[1], D, lb=0, name="zeta_s")
s = SEVRP.addVars(Z_set, D, vertex_ite, lb=0, name="s")

# slack_var = SEVRP.addVars(vertex_ite, D, name="slack_var")

# 'SOCP-related constraints, iteratively enlarged with the number of vertices'
coupling_exprs = [[quicksum(lambda_ce_val[v][i, d] * h_e[d][i] foriinrange(num_eq_cons_ot)) +quicksum(lambda_ci_val[v][j, k, d] * (h_1[d][j] +quicksum(H_3[d][j, i] *Z[i] foriinH_3[d][j].indices) -quicksum(W_b[d][j, i] *x_kad[i, d, k] foriinW_b[d][j].indices)) forjinrange(num_ineq_cons - num_nodes-2-num_cc) forkinK) forvinvertex_ite] fordinD]
Z11_expr = [[gamma[d] -quicksum(zeta_s[i, d]*b[i, d] foriinrange(A.shape[1])) + coupling_exprs[d][v] fordinD] forvinvertex_ite]
Z21_expr = [[[0.5 * (alpha[j, d] +quicksum(A_T[d][j, i]*zeta_s[i, d] foriinrange(A.shape[1])) -quicksum(lambda_ci_val[v][i, k, d]*H_2[d][i, j] foriin (H_2[d].T)[j].indices forkinK)) forjinrange(H_2[d].shape[1])] fordinD] forvinvertex_ite]
Z22 = Q
SEVRP.addConstrs((s.sum('*', d, v) ==Z11_expr[v][d] forvinvertex_itefordinD), name=f'socp_first_cons')
SEVRP.addConstrs((s[j, d, v] +Z22[j, d] >=0forjinZ_setforvinvertex_itefordinD), name=f'socp_trace')

start_time = time.time()
# SEVRP.addConstrs(((s[j, d, v] - Z22[j, d])**2 + 4 * (Z21_expr[v][d][j]) ** 2 <= (s[j, d, v] + Z22[j, d])**2 for j in Z_set for v in vertex_ite for d in D), name=f'socp_det')
SEVRP.addConstrs((s[j, d, v] *Z22[j, d] -Z21_expr[v][d][j]**2 >= 0forjinZ_setforvinvertex_itefordinD), name=f'socp_det')
# for j in Z_set:
# for v in vertex_ite:
# for d in D:
# SEVRP.addQConstr(s[j, d, v] * Z22[j, d] - Z21_expr[v][d][j]**2 >= 0)
end_time = time.time()
print(f"Time taken to add SOCP constraint: {end_time - start_time} seconds")