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.

def eld_complete(U,p_min,p_max,d,brk):
    """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.
    """

    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