mult_selection -- add piecewise relation with multiple selection formulation
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.

def mult_selection(model,a,b):
    """mult_selection -- add piecewise relation with multiple selection formulation
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.
    """

    K = len(a)-1
    w,z = {},{}
    for k in range(K):
        w[k] = model.addVar(lb=-model.infinity()) # do not name variables for avoiding clash
        z[k] = model.addVar(vtype="B")
    X = model.addVar(lb=a[0], ub=a[K], vtype="C")
    Y = model.addVar(lb=-model.infinity())

    for k in range(K):
        model.addCons(w[k] >= a[k]*z[k])
        model.addCons(w[k] <= a[k+1]*z[k])

    model.addCons(quicksum(z[k] for k in range(K)) == 1)
    model.addCons(X == quicksum(w[k] for k in range(K)))

    c = [float(b[k+1]-b[k])/(a[k+1]-a[k]) for k in range(K)]
    d = [b[k]-c[k]*a[k] for k in range(K)]
    model.addCons(Y == quicksum(d[k]*z[k] + c[k]*w[k] for k in range(K)))

    return X,Y,z