code: def eld_complete(U,p_min,p_max,d,brk):
      
  
      model = Model("Economic load dispatching")
  
      p,F = {},{}
      for u in U:
          p[u] = model.addVar(lb=p_min[u], ub=p_max[u], name="p(%s)"%u)    # capacity
          F[u] = model.addVar(lb=0,name="fuel(%s)"%u)
  
      # set fuel costs based on piecewise linear approximation
      for u in U:
          abrk = [X for (X,Y) in brk[u]]
          bbrk = [Y for (X,Y) in brk[u]]
  
          # convex combination part:
          K = len(brk[u])-1
          z = {}
          for k in range(K+1):
              z[k] = model.addVar(ub=1) # do not name variables for avoiding clash
  
          model.addCons(p[u] == quicksum(abrk[k]*z[k] for k in range(K+1)))
          model.addCons(F[u] == quicksum(bbrk[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)])
  
      # demand satisfaction
      model.addCons(quicksum(p[u] for u in U) == d, "demand")
  
      # objective
      model.setObjective(quicksum(F[u] for u in U), "minimize")
  
      model.data = p
      return model
docs: eld -- economic load dispatching in electricity generation
      Parameters:
          - U: set of generators (units)
          - p_min[u]: minimum operating power for unit u
          - p_max[u]: maximum operating power for unit u
          - d: demand
          - brk[k]: (x,y) coordinates of breakpoint k, k=0,...,K
      Returns a model, ready to be solved.