Issue with Mathematical Model Using PuLP Library

I've created a Python script using the PuLP library to solve the multi-skill resource-constrained project scheduling problem (MS-RCPSP), but I'm encountering issues with the task assignments. The output I'm getting doesn't match the expected results shown in Fig-1. Can you help me understand what's going wrong?

import pulp as pl

# Problem Setup
model = pl.LpProblem("Task_Scheduling", pl.LpMinimize)

# Sets
Nr = range(5)  # Tasks
Rr = range(4)  # Resources
Sr = range(2)  # Skills
Kr = range(4)  # Resources skill levels
Ir = range(5)  # Tasks
Jr = range(2)  # Task skills
M = 20

# Data
Prec = [
    [0, 0, 1, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0]
]

p = [5, 10, 5, 10, 5]  # Duration of each task

R = [
    [1, 2],
    [1, 0],
    [1, 1],
    [0, 1],
    [1, 1]
]

B = [
    [1, 0, 1, 1],
    [1, 1, 0, 0]
]

# Decision Variables
x = pl.LpVariable.dicts("x", (Ir, Jr, Kr), cat=pl.LpBinary)
z = pl.LpVariable.dicts("z", (Ir, Ir), cat=pl.LpBinary)
s = pl.LpVariable.dicts("s", Ir, cat=pl.LpContinuous, lowBound=0)
C_max = pl.LpVariable("C_max", lowBound=0, cat=pl.LpContinuous)

# Objective Function
model += C_max, "Minimize_Max_Completion_Time"

# Constraints for max completion time
for j in Nr:
    model += C_max >= s[j] + p[j], f"Max_Completion_Time_{j}"

# Skill requirements
for i in Nr:
    for j in Sr:
        model += pl.lpSum(x[i][j][k] for k in Rr) == R[i][j], f"Skill_Requirement_{i}_{j}"

# Precedence and non-overlap constraints
for i in Nr:
    for i_prime in Nr:
        if i < i_prime:
            model += s[i] + p[i] - M * (1 - z[i][i_prime]) * (1 - Prec[i][i_prime]) <= s[i_prime], f"Precedence_{i}_{i_prime}"
            model += z[i][i_prime] + z[i_prime][i] <= 1, f"Non_Overlap_{i}_{i_prime}"

# Resource constraints
for i in Nr:
    for k in Rr:
        model += pl.lpSum(x[i][j][k] for j in Sr) <= 1, f"One_Skill_Per_Resource_{i}_{k}"

# Resource sharing constraints
for i in Nr:
    for i_prime in Nr:
        if i < i_prime:
            for k in Rr:
                model += pl.lpSum(x[i][j][k] for j in Sr) + pl.lpSum(x[i_prime][j][k] for j in Sr) <= 1 + z[i][i_prime] + z[i_prime][i], f"No_Simultaneous_Assignment_{i}_{i_prime}_{k}"

# Skill level matching
for i in Nr:
    for j in Sr:
        for k in Rr:
            model += x[i][j][k] <= B[j][k], f"Skill_Level_Match_{i}_{j}_{k}"

# Solve the model
model.solve()
print("Status:", pl.LpStatus[model.status])
print("Maximum completion time:", pl.value(C_max))

# Display assignment results and timing for each task
for i in Nr:
    start_time = pl.value(s[i])
    finish_time = start_time + p[i]
    print(f"Activity {i+1} starts at time {start_time} and finishes at time {finish_time}.")
    for j in Sr:
        for k in Rr:
            if pl.value(x[i][j][k]) == 1:
                print(f"Resource {k+1} is assigned to skill {j+1} for activity {i+1}")

I don't think you coded resource requirement constraint (equation 2) correctly.

As I read the code, it does not enforce that resource k has skill j, and thus you are getting an illegal assignment in your results for Activity 4.

So there are 2 ways to do this....

One: multiply the variable by the binary parameter from B like this:

# Skill requirements
for i in Nr:
    for j in Sr:
        model += pl.lpSum(x[i][j][k] * B[j][k] for k in Rr) == R[i][j], f"Skill_Requirement_{i}_{j}"

This is linear, and by inspection, can only credit the resource requirement if B[j][k] is 1.

A second approach would be to make subsets of your x variable with different skills and draw from the appropriate subset to make the equality constraint. Or do this within the summation with a conditional so that only appropriate x are considered. Net effect is the same. Making the subsets might produce a slightly more compact model.

When I did the above, I noticed there were overlapping assignments for resources... I note that in your precedence constraints, I think you shorted the iteration for the main part, equation 3, which should be done for all i, i'. Move the conditional down so it just works on the second equation:

# Precedence and non-overlap constraints
for i in Nr:
    for i_prime in Nr:
        model += s[i] + p[i] - M * (1 - z[i][i_prime]) * (1 - Prec[i][i_prime]) <= s[i_prime], f"Precedence_{i}_{i_prime}"
        if i < i_prime:
            model += z[i][i_prime] + z[i_prime][i] <= 1, f"Non_Overlap_{i}_{i_prime}"

This appears to work and produce a quality answer. It appears there are multiple optimal solutions (as there must be because resource 3 and 4 are interchangeable), as this differs slightly, but appears legit. Note I'm also printing off the z matrix for verification with this addendum to your code:

# Display assignment results and timing for each task
for i in Nr:
    start_time = pl.value(s[i])
    finish_time = start_time + p[i]
    print(f"Activity {i+1} starts at time {start_time} and finishes at time {finish_time}.")
    for j in Sr:
        for k in Rr:
            if pl.value(x[i][j][k]) == 1:
                print(f"Resource {k+1} is assigned to skill {j+1} for activity {i+1}")
    for ii in Nr:
        print(pl.value(z[i][ii]), end = ' ')
    print()
    print()

Output:

Status: Optimal
Maximum completion time: 15.0
Activity 1 starts at time 0.0 and finishes at time 5.0.
Resource 4 is assigned to skill 1 for activity 1
Resource 1 is assigned to skill 2 for activity 1
Resource 2 is assigned to skill 2 for activity 1
0.0 0.0 1.0 1.0 1.0 

Activity 2 starts at time 0.0 and finishes at time 10.0.
Resource 3 is assigned to skill 1 for activity 2
0.0 0.0 0.0 0.0 1.0 

Activity 3 starts at time 5.0 and finishes at time 10.0.
Resource 4 is assigned to skill 1 for activity 3
Resource 1 is assigned to skill 2 for activity 3
0.0 0.0 0.0 0.0 1.0 

Activity 4 starts at time 5.0 and finishes at time 15.0.
Resource 2 is assigned to skill 2 for activity 4
0.0 0.0 0.0 0.0 0.0 

Activity 5 starts at time 10.0 and finishes at time 15.0.
Resource 4 is assigned to skill 1 for activity 5
Resource 1 is assigned to skill 2 for activity 5
0.0 0.0 0.0 0.0 0.0