rcs -- model for the resource constrained scheduling problem
    Parameters:
        - J: set of jobs
        - P: set of precedence constraints between jobs
        - R: set of resources
        - T: number of periods
        - p[j]: processing time of job j
        - c[j,t]: cost incurred when job j starts processing on period t.
        - a[j,r,t]: resource r usage for job j on period t (after job starts)
        - RUB[r,t]: upper bound for resource r on period t
    Returns a model, ready to be solved.

def rcs(J,P,R,T,p,c,a,RUB):
    """rcs -- model for the resource constrained scheduling problem
    Parameters:
        - J: set of jobs
        - P: set of precedence constraints between jobs
        - R: set of resources
        - T: number of periods
        - p[j]: processing time of job j
        - c[j,t]: cost incurred when job j starts processing on period t.
        - a[j,r,t]: resource r usage for job j on period t (after job starts)
        - RUB[r,t]: upper bound for resource r on period t
    Returns a model, ready to be solved.
    """
    model = Model("resource constrained scheduling")

    s,x = {},{}   # s - start time variable;  x=1 if job j starts on period t
    for j in J:
        s[j] = model.addVar(vtype="C", name="s(%s)"%j)
        for t in range(1,T-p[j]+2):
            x[j,t] = model.addVar(vtype="B", name="x(%s,%s)"%(j,t))

    for j in J:
        # job execution constraints
        model.addCons(quicksum(x[j,t] for t in range(1,T-p[j]+2)) == 1, "ConstrJob(%s,%s)"%(j,t))

        # start time constraints
        model.addCons(quicksum((t-1)*x[j,t] for t in range(2,T-p[j]+2)) == s[j], "ConstrJob(%s,%s)"%(j,t))

    # resource upper bound constraints
    for t in range(1,T-p[j]+2):
        for r in R:
            model.addCons(
                quicksum(a[j,r,t-t_]*x[j,t_] for j in J for t_ in range(max(t-p[j]+1,1),min(t+1,T-p[j]+2))) \
                <= RUB[r,t], "ResourceUB(%s)"%t)

    # time (precedence) constraints, i.e., s[k]-s[j] >= p[j]
    for (j,k) in P:
        model.addCons(s[k] - s[j] >= p[j], "Precedence(%s,%s)"%(j,k))

    model.setObjective(quicksum(c[j,t]*x[j,t] for (j,t) in x), "minimize")

    model.data = x,s
    return model