code: def mult_selection(model,a,b):
      
  
      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
docs: 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.