Problem with solution of soc relaxation of an optimal power flow problem

Hello,

I implemented (hardcoded) the following jabr-soc-relaxation formulation of the opf from https://ieeexplore.ieee.org/document/7056568 using 2-bus-2-generator instances from IV A:

import numpy as np
import gurobipy as gp
from gurobipy import GRB


def sdp_power_flow_optimization_gurobi(
    num_buses,
    generator_index,
    costs,
    L,
    p_d,
    q_d,
    G_ij,
    B_ij,
    V_min,
    V_max,
    p_min,
    p_max,
    q_min,
    q_max,
):
    # Create a new model
    m = gp.Model("power_flow")

    # Define variables
    c_ij = {}
    s_ij = {}
    # c_ii = {}
    for i, j in L:
        c_ij[(i, j)] = m.addVar(name=f"c_{i}{j}")
        s_ij[(i, j)] = m.addVar(name=f"s_{i}{j}")

    """
    for i in range(num_buses):
        c_ii[i] = m.addVar(name=f"c_{i}{i}")
    """
    p_g = m.addVars(num_buses, name="p_g")
    q_g = m.addVars(num_buses, name="q_g")

    base_unit = 100

    # Set objective, linear generator cost func
    m.setObjective(
        base_unit * gp.quicksum(costs[i] * p_g[i] for i in range(num_buses)),
        GRB.MINIMIZE,
    )

    # Add constraints
    # Power balance constraints

    # Set hardcoded shunt elements and voltage magnitudes from MATPOWER
    G_12 = -3.8156
    B_12 = 19.0781
    G_11 = 3.8156
    B_11 = -19.0781

    G_21 = -3.8156
    B_21 = 19.0781
    G_22 = 3.8156
    B_22 = -19.0781

    c_11 = 0.874
    c_22 = 0.816

    slack_p = m.addVars(num_buses, name="slack_p", lb=-GRB.INFINITY)
    slack_q = m.addVars(num_buses, name="slack_q", lb=-GRB.INFINITY)

    # Hardcoded power-balance constraints
    m.addConstr(
        p_g[0] - (75 / base_unit) + slack_p[0]
        == G_12 * c_ij[(0, 1)] - B_12 * s_ij[(0, 1)] + G_11 * c_11
    )
    m.addConstr(
        q_g[0] + (84.7 / base_unit) + slack_q[0]
        == -B_12 * c_ij[(0, 1)] - G_12 * s_ij[(0, 1)] - B_11 * c_11
    )
    m.addConstr(
        p_g[1] - (105 / base_unit) + slack_p[1]
        == G_21 * c_ij[(1, 0)] - B_21 * s_ij[(1, 0)] + G_22 * c_22
    )
    m.addConstr(
        q_g[1] - (22.8 / base_unit) + slack_q[1]
        == -B_21 * c_ij[(1, 0)] - G_21 * s_ij[(1, 0)] - B_22 * c_22
    )

    """
    # Voltage magnitude constraints
    for i in range(num_buses):
        m.addConstr(c_ii[i] >= V_min[i] ** 2)
        m.addConstr(c_ii[i] <= V_max[i] ** 2)
    """

    # Generation limits
    m.addConstr(p_g[0] >= p_min[0])
    m.addConstr(p_g[0] <= p_max[0])
    m.addConstr(q_g[0] >= q_min[0])
    m.addConstr(q_g[0] <= q_max[0])

    m.addConstr(p_g[1] >= p_min[1])
    m.addConstr(p_g[1] <= p_max[1])
    m.addConstr(q_g[1] >= q_min[1])
    m.addConstr(q_g[1] <= q_max[1])

    m.addConstr(c_ij[(0, 1)] == c_ij[(1, 0)])
    m.addConstr(s_ij[(0, 1)] == -s_ij[(1, 0)])

    # Line constraints
    m.addQConstr(c_ij[(0, 1)] ** 2 + s_ij[(0, 1)] ** 2 <= c_11 * c_22)

    # Optimize the model
    m.setParam("NonConvex", 2)
    m.optimize()

    # Check if the model is infeasible
    if m.status == GRB.INF_OR_UNBD or m.status == GRB.INFEASIBLE:
        print("Model is infeasible")

        # Compute the IIS
        m.computeIIS()

        # Print the IIS
        print("\nThe following constraint(s) cannot be satisfied:")
        for c in m.getConstrs():
            if c.IISConstr:
                # Get the constraint expression (LHS)
                lhs_expr = m.getRow(c)

                # Get the RHS of the constraint
                sense = c.Sense
                rhs = c.RHS

                # Reconstruct the full constraint expression
                constraint_str = f"{lhs_expr} {sense} {rhs}"
                print(f"{c.constrName}: {constraint_str}")

        # For variables, you can check for IISLB and IISUB attributes
        for v in m.getVars():
            if v.IISLB > 0 or v.IISUB > 0:
                print("Variable bounds causing issue: %s" % v.varName)

    if m.status == GRB.OPTIMAL:
        # Iterate over the variables and print their names and optimized values
        for v in m.getVars():
            print(f"{v.varName} = {v.x}")
    else:
        print("Optimal solution was not found.")

    return 0


# Example from Inexactness of SDP Relaxation ... Kocuk, Dey and Sun
num_buses = 2
generator_index = {0, 1}

costs = np.array([5, 1.2])  # Here costs per MW
p_d = np.array([75, 105])  # Example demand power
q_d = np.array([-84.7, 22.8])  # Example demand reactive power

# Set of transmission lines
L = {(0, 1), (1, 0)}
# shunt conductance
g_ii = np.array([0, 0])
# shunt susceptance
b_ii = np.array([0, 0])

base_unit = 100

V_min = np.array([90 / base_unit, 90 / base_unit])
V_max = np.array([110 / base_unit, 110 / base_unit])
p_min = np.array([75 / base_unit, 70 / base_unit])
p_max = np.array([250 / base_unit, 300 / base_unit])
q_min = np.array([-30 / base_unit, -30 / base_unit])
q_max = np.array([300 / base_unit, 300 / base_unit])

# Call the optimization function
results = sdp_power_flow_optimization_gurobi(
    num_buses,
    generator_index,
    costs,
    L,
    p_d,
    q_d,
    G_ij,
    B_ij,
    V_min,
    V_max,
    p_min,
    p_max,
    q_min,
    q_max,
)

Now I am facing the issue that the optimal objective values are not the same and if I change the soc constraint to:

m.addQConstr(c_ij[(0, 1)] ** 2 + s_ij[(0, 1)] ** 2 == c_11 * c_22)

the problem becomes infeasible which should not happen. I am not sure what to do or how to adjust the problem so that I can generate same results. 
I was thinking that probably gurobi generates different solutions than for instance mosek but I think problems of infeasibility doesnt depend on the solver for such easy instances. 

Very, very grateful for any hints/support

Paul