code: def convex_comb_sos(model,a,b):
      
      K = len(a)-1
      z = {}
      for k in range(K+1):
          z[k] = model.addVar(lb=0, ub=1, vtype="C")
      X = model.addVar(lb=a[0], ub=a[K], vtype="C")
      Y = model.addVar(lb=-model.infinity(), vtype="C")
  
      model.addCons(X == quicksum(a[k]*z[k] for k in range(K+1)))
      model.addCons(Y == quicksum(b[k]*z[k] for k in range(K+1)))
  
      model.addCons(quicksum(z[k] for k in range(K+1)) == 1)
      model.addConsSOS2([z[k] for k in range(K+1)])
  
      return X,Y,z
docs: convex_comb_sos -- add piecewise relation with gurobi's SOS constraints
      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.