weber -- model for solving the weber problem using soco (multiple source version).
    Parameters:
        - I: set of customers
        - J: set of potential facilities
        - x[i]: x position of customer i
        - y[i]: y position of customer i
        - w[i]: weight of customer i
    Returns a model, ready to be solved.

def weber_MS(I,J,x,y,w):
    """weber -- model for solving the weber problem using soco (multiple source version).
    Parameters:
        - I: set of customers
        - J: set of potential facilities
        - x[i]: x position of customer i
        - y[i]: y position of customer i
        - w[i]: weight of customer i
    Returns a model, ready to be solved.
    """
    M = max([((x[i]-x[j])**2 + (y[i]-y[j])**2) for i in I for j in I])
    model = Model("weber - multiple source")
    X,Y,v,u = {},{},{},{}
    xaux,yaux,uaux = {},{},{}
    for j in J:
        X[j] = model.addVar(lb=-model.infinity(), vtype="C", name="X(%s)"%j)
        Y[j] = model.addVar(lb=-model.infinity(), vtype="C", name="Y(%s)"%j)
        for i in I:
            v[i,j] = model.addVar(vtype="C", name="v(%s,%s)"%(i,j))
            u[i,j] = model.addVar(vtype="B", name="u(%s,%s)"%(i,j))
            xaux[i,j] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s,%s)"%(i,j))
            yaux[i,j] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s,%s)"%(i,j))
            uaux[i,j] = model.addVar(vtype="C", name="uaux(%s,%s)"%(i,j))



    for i in I:
        model.addCons(quicksum(u[i,j] for j in J) == 1, "Assign(%s)"%i)
        for j in J:
            model.addCons(xaux[i,j]*xaux[i,j] + yaux[i,j]*yaux[i,j] <= v[i,j]*v[i,j], "MinDist(%s,%s)"%(i,j))
            model.addCons(xaux[i,j] == (x[i]-X[j]), "xAux(%s,%s)"%(i,j))
            model.addCons(yaux[i,j] == (y[i]-Y[j]), "yAux(%s,%s)"%(i,j))
            model.addCons(uaux[i,j] >= v[i,j] - M*(1-u[i,j]), "uAux(%s,%s)"%(i,j))

    model.setObjective(quicksum(w[i]*uaux[i,j] for i in I for j in J), "minimize")


    model.data = X,Y,v,u
    return model