Trying to understand why my formulation is non-convex? (gurobipy.GurobiError: Q matrix is not positive semi-definite (PSD))

Question

I am trying to implement a SOCP that I believe to be convex, but when I run it on a small instance I get the error:

gurobipy.GurobiError: Q matrix is not positive semi-definite (PSD). Set NonConvex parameter to 2 to solve model.

Does anyone have any advice on how to understand why it is non-convex?

Here is the small instance of the SOCP model for 3 vertices, 6 edges:

Objective: 0.81*chi_m['1'] - 1.2100000000000002*chi_p['1'] - 1.1*(tau['1'] + xi['1']) - 0.4*(eta['1'] + zeta['1']) + (0.81*chi_m['2'] - 1.2100000000000002*chi_p['2']) - 1.1*(tau['2'] + xi['2']) - 0.4*(eta['2'] + zeta['2']) + (0.81*chi_m['3'] - 1.2100000000000002*chi_p['3']) - 0.95*(tau['3'] + xi['3']) - 0.5*(eta['3'] + zeta['3']) + (0.0*gamma_m['1'] - 20.0*gamma_p['1'] - 10.0*kappa_m['1'] - 10.0*kappa_p['1'] + 0.0*gamma_m['2'] - 20.0*gamma_p['2'] - 10.0*kappa_m['2'] - 10.0*kappa_p['2'] + 0.0*gamma_m['3'] - 0.0*gamma_p['3'] - 10.0*kappa_m['3'] - 10.0*kappa_p['3']) - (8100.0*omega['1'] + 0.25*omega['2'] + 8100.0*omega['3'] + 8100.0*omega['4'] + 0.25*omega['5'] + 8100.0*omega['6'])

power_flow_conic_dual_constr['1']: epsilon['1',0]**2 + epsilon['1',1]**2  <=  omega['1']**2
power_flow_conic_dual_constr['2']: epsilon['2',0]**2 + epsilon['2',1]**2  <=  omega['2']**2
power_flow_conic_dual_constr['3']: epsilon['3',0]**2 + epsilon['3',1]**2  <=  omega['3']**2
power_flow_conic_dual_constr['4']: epsilon['4',0]**2 + epsilon['4',1]**2  <=  omega['4']**2
power_flow_conic_dual_constr['5']: epsilon['5',0]**2 + epsilon['5',1]**2  <=  omega['5']**2
power_flow_conic_dual_constr['6']: epsilon['6',0]**2 + epsilon['6',1]**2  <=  omega['6']**2

relaxation_conic_dual_constr['1']: delta['1',0]**2 + delta['1',1]**2 + delta['1',2]**2 + delta['1',3]**2  <=  sigma['1']**2
relaxation_conic_dual_constr['2']: delta['2',0]**2 + delta['2',1]**2 + delta['2',2]**2 + delta['2',3]**2  <=  sigma['2']**2
relaxation_conic_dual_constr['3']: delta['3',0]**2 + delta['3',1]**2 + delta['3',2]**2 + delta['3',3]**2  <=  sigma['3']**2
relaxation_conic_dual_constr['4']: delta['4',0]**2 + delta['4',1]**2 + delta['4',2]**2 + delta['4',3]**2  <=  sigma['4']**2
relaxation_conic_dual_constr['5']: delta['5',0]**2 + delta['5',1]**2 + delta['5',2]**2 + delta['5',3]**2  <=  sigma['5']**2
relaxation_conic_dual_constr['6']: delta['6',0]**2 + delta['6',1]**2 + delta['6',2]**2 + delta['6',3]**2  <=  sigma['6']**2

lift_e_dual_constr['1']: chi_m['1'] - chi_p['1'] + 0.0*xi['1'] + 0.0*zeta['1'] - (0.16725635252492763*phi['1'] + 1.3703682856223864*psi['1'] + sigma['1'] + delta['1',2] + 0.05173917542536994*phi['3'] + 0.9586966162579272*psi['3'] + sigma['3'] + delta['3',2]) - (sigma['4'] - delta['4',3] + sigma['6'] - delta['6',3])  <=  0

lift_e_dual_constr['2']: chi_m['2'] - chi_p['2'] + 0.0*xi['2'] + 0.0*zeta['2'] - (0.04439511653718091*phi['5'] + 0.9818534961154274*psi['5'] + sigma['5'] + delta['5',2] + 0.05173917542536994*phi['6'] + 0.9586966162579272*psi['6'] + sigma['6'] + delta['6',2]) - (sigma['2'] - delta['2',3] + sigma['3'] - delta['3',3])  <=  0

lift_e_dual_constr['3']: chi_m['3'] - chi_p['3'] + 0.0*xi['3'] + 0.0*zeta['3'] - (0.04439511653718091*phi['2'] + 0.9818534961154274*psi['2'] + sigma['2'] + delta['2',2] + 0.16725635252492763*phi['4'] + 1.3703682856223864*psi['4'] + sigma['4'] + delta['4',2]) - (sigma['1'] - delta['1',3] + sigma['5'] - delta['5',3])  <=  0

pg_dual_constr['1']: gamma_m['1'] - gamma_p['1'] + xi['1']  ==  0
pg_dual_constr['2']: gamma_m['2'] - gamma_p['2'] + xi['2']  ==  0
pg_dual_constr['3']: gamma_m['3'] - gamma_p['3'] + xi['3']  ==  0

qg_dual_constr['1']: kappa_m['1'] - kappa_p['1'] + zeta['1']  ==  0
qg_dual_constr['2']: kappa_m['2'] - kappa_p['2'] + zeta['2']  ==  0
qg_dual_constr['3']: kappa_m['3'] - kappa_p['3'] + zeta['3']  ==  0

pu_dual_constr['1']: xi['1'] - tau['1']  <=  1
pu_dual_constr['2']: xi['2'] - tau['2']  <=  1
pu_dual_constr['3']: xi['3'] - tau['3']  <=  1

qu_dual_constr['1']: zeta['1'] - eta['1']  <=  0
qu_dual_constr['2']: zeta['2'] - eta['2']  <=  0
qu_dual_constr['3']: zeta['3'] - eta['3']  <=  0

pij_dual_constr['1']: xi['1'] + xi['3'] + phi['1'] - epsilon['1',0]  ==  0
pij_dual_constr['2']: xi['3'] + xi['2'] + phi['2'] - epsilon['2',0]  ==  0
pij_dual_constr['3']: xi['1'] + xi['2'] + phi['3'] - epsilon['3',0]  ==  0
pij_dual_constr['4']: xi['3'] + xi['1'] + phi['4'] - epsilon['4',0]  ==  0
pij_dual_constr['5']: xi['2'] + xi['3'] + phi['5'] - epsilon['5',0]  ==  0
pij_dual_constr['6']: xi['2'] + xi['1'] + phi['6'] - epsilon['6',0]  ==  0

qij_dual_constr['1']: zeta['1'] + zeta['3'] + psi['1'] - epsilon['1',1]  ==  0
qij_dual_constr['2']: zeta['3'] + zeta['2'] + psi['2'] - epsilon['2',1]  ==  0
qij_dual_constr['3']: zeta['1'] + zeta['2'] + psi['3'] - epsilon['3',1]  ==  0
qij_dual_constr['4']: zeta['3'] + zeta['1'] + psi['4'] - epsilon['4',1]  ==  0
qij_dual_constr['5']: zeta['2'] + zeta['3'] + psi['5'] - epsilon['5',1]  ==  0
qij_dual_constr['6']: zeta['2'] + zeta['1'] + psi['6'] - epsilon['6',1]  ==  0

lift_s_dual_constr['1']: -1.5953682856223865*phi['1'] + 0.16725635252492763*psi['1'] - 2*delta['1',0] + pi['1'] + pi['4']  ==  0     
lift_s_dual_constr['2']: -1.3318534961154274*phi['2'] + 0.04439511653718091*psi['2'] - 2*delta['2',0] + pi['2'] + pi['5']  ==  0     
lift_s_dual_constr['3']: -1.1086966162579273*phi['3'] + 0.05173917542536994*psi['3'] - 2*delta['3',0] + pi['3'] + pi['6']  ==  0     
lift_s_dual_constr['4']: -1.5953682856223865*phi['4'] + 0.16725635252492763*psi['4'] - 2*delta['4',0] + pi['4'] + pi['1']  ==  0     
lift_s_dual_constr['5']: -1.3318534961154274*phi['5'] + 0.04439511653718091*psi['5'] - 2*delta['5',0] + pi['5'] + pi['2']  ==  0     
lift_s_dual_constr['6']: -1.1086966162579273*phi['6'] + 0.05173917542536994*psi['6'] - 2*delta['6',0] + pi['6'] + pi['3']  ==  0     

lift_c_dual_constr['1']: 0.16725635252492763*phi['1'] + 1.5953682856223865*psi['1'] - 2*delta['1',1] + rho['1'] - rho['4']  ==  0    
lift_c_dual_constr['2']: 0.04439511653718091*phi['2'] + 1.3318534961154274*psi['2'] - 2*delta['2',1] + rho['2'] - rho['5']  ==  0    
lift_c_dual_constr['3']: 0.05173917542536994*phi['3'] + 1.1086966162579273*psi['3'] - 2*delta['3',1] + rho['3'] - rho['6']  ==  0    
lift_c_dual_constr['4']: 0.16725635252492763*phi['4'] + 1.5953682856223865*psi['4'] - 2*delta['4',1] + rho['4'] - rho['1']  ==  0    
lift_c_dual_constr['5']: 0.04439511653718091*phi['5'] + 1.3318534961154274*psi['5'] - 2*delta['5',1] + rho['5'] - rho['2']  ==  0    
lift_c_dual_constr['6']: 0.05173917542536994*phi['6'] + 1.1086966162579273*psi['6'] - 2*delta['6',1] + rho['6'] - rho['3']  ==  0    

The variable bounds are as follows:

model.gamma_p = Var(gen_list, within=NonNegativeReals) #  γ^+ 
model.gamma_m = Var(gen_list, within=NonNegativeReals) #  γ^-
model.kappa_p = Var(gen_list, within=NonNegativeReals) # κ^+
model.kappa_m = Var(gen_list, within=NonNegativeReals) # κ^-
model.chi_p = Var(bus_list, within=NonNegativeReals) # χ^+
model.chi_m = Var(bus_list, within=NonNegativeReals) # χ^-
model.tau = Var(bus_list, within=NonNegativeReals) # τ
model.eta = Var(bus_list, within=NonNegativeReals) # η
model.xi = Var(bus_list, within=Reals) # ξ
model.zeta = Var(bus_list, within=Reals) # ζ
model.phi = Var(branch_list, within=Reals) # φ
model.psi = Var(branch_list, within=Reals) # ψ
model.omega = Var(branch_list, within=Reals) # ω
model.delta = Var(branch_list,[0,1,2], within=Reals) # δ
model.epsilon = Var(branch_list,[0,1], within=Reals) # ε
model.sigma = Var(branch_list, within=Reals) # σ
model.pi = Var(branch_list, within=Reals) # π
model.rho = Var(branch_list, within=Reals) # ρ

Answer 1

The input format of a second-order cone constraint looks like:

x^2 + y^2 <= z^2
with z a non-negative variable

You are not following this. E.g. you have

epsilon['1',0]**2 + epsilon['1',1]**2  <=  omega['1']**2

with all free variables.

This is actually properly documented: see e.g. https://www.gurobi.com/documentation/10.0/refman/constraints.html.