ssa -- multi-stage (serial) safety stock allocation model
    Parameters:
        - n: number of stages
        - h[i]: inventory cost on stage i
        - K: number of linear segments
        - f: (non-linear) cost function
        - T[i]: production lead time on stage i
    Returns the model with the piecewise linear relation on added variables x, f, and z.

def ssa(n,h,K,f,T):
    """ssa -- multi-stage (serial) safety stock allocation model
    Parameters:
        - n: number of stages
        - h[i]: inventory cost on stage i
        - K: number of linear segments
        - f: (non-linear) cost function
        - T[i]: production lead time on stage i
    Returns the model with the piecewise linear relation on added variables x, f, and z.
    """

    model = Model("safety stock allocation")

    # calculate endpoints for linear segments
    a,b = {},{}
    for i in range(1,n+1):
        a[i] = [k for k in range(K)]
        b[i] = [f(i,k) for k in range(K)]

    # x: net replenishment time for stage i
    # y: corresponding cost
    # s: piecewise linear segment of variable x
    x,y,s = {},{},{}
    L = {} # service time of stage i
    for i in range(1,n+1):
        x[i],y[i],s[i] = convex_comb_sos(model,a[i],b[i])
        if i == 1:
            L[i] = model.addVar(ub=0, vtype="C", name="L[%s]"%i)
        else:
            L[i] = model.addVar(vtype="C", name="L[%s]"%i)
    L[n+1] = model.addVar(ub=0, vtype="C", name="L[%s]"%(n+1))

    for i in range(1,n+1):
        # net replenishment time for each stage i
        model.addCons(x[i] + L[i] == T[i] + L[i+1])

    model.setObjective(quicksum(h[i]*y[i] for i in range(1,n+1)), "minimize")

    model.data = x,s,L
    return model