flp_nonlinear_mselect --  use multiple selection model
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for customer 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
        - K: number of linear pieces for approximation of non-linear cost function
    Returns a model, ready to be solved.

def flp_nonlinear_mselect(I,J,d,M,f,c,K):
    """flp_nonlinear_mselect --  use multiple selection model
    Parameters:
        - I: set of customers
        - J: set of facilities
        - d[i]: demand for customer 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
        - K: number of linear pieces for approximation of non-linear cost function
    Returns a model, ready to be solved.
    """
    a,b = {},{}
    for j in J:
        U = M[j]
        L = 0
        width = U/float(K)
        a[j] = [k*width for k in range(K+1)]
        b[j] = [f[j]*math.sqrt(value) for value in a[j]]

    model = Model("nonlinear flp -- piecewise linear version with multiple selection")
    
    x = {}
    for j in J:
        for i in I:
            x[i,j] = model.addVar(vtype="C", name="x(%s,%s)"%(i,j))  # i's demand satisfied from j

    # total volume transported from plant j, corresponding (linearized) cost, selection variable:
    X,F,z = {},{},{}
    for j in J:
        # add constraints for linking piecewise linear part:
        X[j],F[j],z[j] = mult_selection(model,a[j],b[j])
        X[j].ub = M[j]
        # for i in I:
        #     model.addCons(
        #         x[i,j] <= \
        #         quicksum(min(d[i],a[j][k+1]) * z[j][k] for k in range(K)),\
        #         "Strong(%s,%s)"%(i,j))

    # constraints for customer's demand satisfaction
    for i in I:
        model.addCons(quicksum(x[i,j] for j in J) == d[i], "Demand(%s)"%i)

    for j in J:
        model.addCons(quicksum(x[i,j] for i in I) == X[j], "Capacity(%s)"%j)

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

    model.data = x,X,F
    return model