flp_nonlinear_soco --  use
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for product i
        - M[j]: capacity of facility j
        - f[j]: fixed cost for using a facility in point j
        - c[i,j]: unit cost of servicing demand point i from facility j
    Returns a model, ready to be solved.

def flp_nonlinear_soco(I,J,d,M,f,c):
    """flp_nonlinear_soco --  use
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for product i
        - M[j]: capacity of facility j
        - f[j]: fixed cost for using a facility in point j
        - c[i,j]: unit cost of servicing demand point i from facility j
    Returns a model, ready to be solved.
    """
    model = Model("nonlinear flp -- soco formulation")

    x,X,u = {},{},{}
    for j in J:
        X[j] = model.addVar(ub=M[j], vtype="C", name="X(%s)"%j) # for sum_i x_ij
        u[j] = model.addVar(vtype="C", name="u(%s)"%j) # for replacing sqrt sum_i x_ij in soco
        for i in I:
            x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))  # i's demand satisfied from j
    
    # constraints for customer's demand satisfaction
    for i in I:
        model.addCons(quicksum(x[i,j] for j in J) == 1, "Demand(%s)"%i)

    for j in J:
        model.addCons(quicksum(d[i]*x[i,j] for i in I) == X[j], "Capacity(%s)"%j)
        model.addQConstr(quicksum(f[j]**2*d[i]*x[i,j]*x[i,j] for i in I) <= u[j]*u[j], "SOC(%s)"%j)

    model.setObjective(quicksum(u[j] for j in J) +\
                       quicksum(c[i,j]*d[i]*x[i,j] for j in J for i in I),\
                       "minimize")

    model.data = x,u
    return model