scheduling_cutting_plane: heuristic to one machine weighted completion time based on cutting planes

    Use a cutting-plane method as a heuristics for solving the
    one-machine total weighted completion time problem with release
    times.

    Parameters:
        - J: set of jobs
        - p[j]: processing time of job j
        - r[j]: earliest start time of job j
        - w[j]: weighted of job j; the objective is the sum of the weighted completion time

    Returns a tuple with values of:
        - bestC: corresponding completion time of the best found solution
        - bestobj: best, best, hobj

def scheduling_cutting_plane(J,p,r,w):
    """
    scheduling_cutting_plane: heuristic to one machine weighted completion time based on cutting planes

    Use a cutting-plane method as a heuristics for solving the
    one-machine total weighted completion time problem with release
    times.

    Parameters:
        - J: set of jobs
        - p[j]: processing time of job j
        - r[j]: earliest start time of job j
        - w[j]: weighted of job j; the objective is the sum of the weighted completion time

    Returns a tuple with values of:
        - bestC: corresponding completion time of the best found solution
        - bestobj: best, best, hobj
    """

    model = Model("scheduling: cutting plane")

    C = {}   # completion time variable
    for j in J:
        C[j] = model.addVar(lb=r[j]+p[j], obj=w[j], vtype="C", name="C(%s)"%j)


    sumP = sum([p[j] for j in J])
    sumP2 = sum([p[j]**2 for j in J])
    model.addCons(C[j]*p[j] >= sumP2*0.5 + (sumP**2)*0.5, "AllJobs")

    model.setObjective(quicksum(w[j]*C[j] for j in J), "minimize")

    cut = 0
    bestobj = float("inf")
    while True:
        model.optimize()
        sol = sorted([(model.getVal(C[j]),j) for j in C])
        seq = [j for (completion,j) in sol]

        # print("Opt.value=",model.getObjVal()
        # print("current solution:", seq
        hC,hobj = evaluate(seq,p,r,w)
        if hobj < bestobj:
            # print("\t*** updating  best found solution:", bestobj, "--->", hobj, "***"
            bestC = hC
            bestobj = hobj
            best = list(seq)

        S,S_ = [],[]
        sumP,sumP2,f = 0,0,0
        for (C_,j) in sol:
            pj = p[j]
            delta = pj**2 + ((sumP+pj)**2 - sumP**2) - 2*pj*C_
            if f > 0 and delta <= 0:
                S = S_
                break
            f += delta/2.
            sumP2 += pj**2
            sumP += pj
            S_.append(j)

        if S == []:
            break

        cut += 1
        model.addCons(quicksum(C[j]*p[j] for j in S) >= sumP2*0.5 + (sumP**2)*0.5, "Cut(%s)"%cut)

    return bestC,bestobj,best