NSGA Function IndexError: index 3 is out of bounds for axis 0 with size 3

[jfortizl]

hi! I am terrible programing and I reached a mistake in my code that I am just not able to solve at all. I am doing an optimization with pymoo using NSGA 2, but when I run the code i have the error of index out of bounds. can anyone help me to solve this problem?

my code is:
import pymoo
import numpy as np
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.operators.crossover.pntx import TwoPointCrossover
from pymoo.operators.mutation.bitflip import BitflipMutation
from pymoo.operators.sampling.rnd import BinaryRandomSampling
from pymoo.optimize import minimize
from pymoo.core.problem import Problem

Problem configuration

NUMBER_OF_JOBS = 3
NUMBER_OF_MACHINES = 3
T = np.array([10, 20, 30]) # Time targets for each job
gamma_p = 1.0 # Some weight
Q_p = 5 # Some constant

Energy consumption parameters

P = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES, NUMBER_OF_MACHINES) # Power consumption
y = np.random.randint(0, 2, (NUMBER_OF_JOBS, NUMBER_OF_MACHINES, NUMBER_OF_MACHINES)) # Machine selection binary variables
lambda_k = 0.1 # Weight for power consumption
W = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES) # Weight of jobs
eta_k = 0.05 # Weight for job weights
R = np.array([100, 150, 200]) # Resource availability
U_k = 50 # Resource usage
theta_k = 0.2 # Weight for resources
F = 10 # Some factor

Constraint related parameters

U = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES)
S = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES)
H = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES)
KT = np.random.rand(NUMBER_OF_JOBS, NUMBER_OF_MACHINES)
M = 1000 # Large number for constraints
z = np.random.randint(0, 2, (NUMBER_OF_JOBS, NUMBER_OF_MACHINES, NUMBER_OF_JOBS, NUMBER_OF_MACHINES, NUMBER_OF_MACHINES)) # More binary variables

class MyProblem(Problem):
def init(self):
super().init(n_var=NUMBER_OF_JOBS * NUMBER_OF_MACHINES,
n_obj=4,
n_constr=8, # Number of constraints
xl=0,
xu=1)

def _evaluate(self, x, out, *args, **kwargs):
    # Extracting job schedules from individual
    job_schedules = [x[i:i+NUMBER_OF_MACHINES] for i in range(0, len(x), NUMBER_OF_MACHINES)]

    # Calculate objectives
    f1 = np.max([np.sum(job) for job in job_schedules])
    f2 = np.sum([np.max(np.sum(job) - T[i], 0) for i, job in enumerate(job_schedules)])
    f3 = np.sum([np.sum(job) + gamma_p * Q_p for job in job_schedules])
    E = np.sum([np.sum([np.sum([P[i][j][k] * y[i][j][k] * lambda_k for k in range(NUMBER_OF_MACHINES)]) for j in range(len(job))]) for i, job in enumerate(job_schedules)])
    G = np.sum([np.sum([np.sum([W[i][k] * y[i][j][k] * eta_k for k in range(NUMBER_OF_MACHINES)]) for j in range(len(job))]) for i, job in enumerate(job_schedules)])
    B = np.sum([(R[k] - U_k - np.sum([np.sum([P[i][j][k] * y[i][j][k] for j in range(len(job))]) for i, job in enumerate(job_schedules)]) - np.sum([np.sum([W[i][k] * y[i][j][k] for j in range(len(job))]) for i, job in enumerate(job_schedules)])) * theta_k for k in range(NUMBER_OF_MACHINES)])
    f4 = f1 * F + E + B + G

    # Constraint values
    constr_values = []

    # Equation 8 constraint
    constr_values.append(np.sum([U[i] - np.sum([W[i][k] * job_schedules[i][j][k] for k in range(NUMBER_OF_MACHINES)]) - R[i][j] for i in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES)]))

    # Equation 9 constraint
    constr_values.append(np.sum([R[i][j] - np.sum([job_schedules[i][j][k] * P[i][j][k] for k in range(NUMBER_OF_MACHINES)]) - S[i][j] for i in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES)]))

    # Equation 10 constraint
    constr_values.extend([S[i][j] - R[i][j-1] for i in range(NUMBER_OF_JOBS) for j in range(1, NUMBER_OF_MACHINES)])

    # Equation 11 constraint
    constr_values.extend([H[i][j] - U[i][j-1] for i in range(NUMBER_OF_JOBS) for j in range(1, NUMBER_OF_MACHINES)])

    # Equation 12 constraint
    constr_values.extend([S[h][g] - S[i][j] + M * (1 - z[i][j][h][g][k]) - np.sum([job_schedules[i][j][k] * P[i][j][k] for k in range(NUMBER_OF_MACHINES)]) for i in range(NUMBER_OF_JOBS) for h in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES) for g in range(NUMBER_OF_MACHINES) for k in range(NUMBER_OF_MACHINES)])

    # Equation 13 constraint
    constr_values.extend([H[h][g] - KT[i][j] + M * (1 - z[i][j][h][g][k]) - np.sum([job_schedules[i][j][k] * P[i][j][k] for k in range(NUMBER_OF_MACHINES)]) for i in range(NUMBER_OF_JOBS) for h in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES) for g in range(NUMBER_OF_MACHINES) for k in range(NUMBER_OF_MACHINES)])

    # Equation 14 constraint
    constr_values.extend([np.sum([job_schedules[i][j][k] for k in range(NUMBER_OF_MACHINES)]) - 1 for i in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES)])

    # Constraint values for binary variables (equation 15)
    constr_values.extend([job_schedules[i][j][k] * (1 - job_schedules[i][j][k]) for i in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES) for k in range(NUMBER_OF_MACHINES)])

    # Constraint values for binary variables (equation 16)
    constr_values.extend([z[i][j][h][g][k] * (1 - z[i][j][h][g][k]) for i in range(NUMBER_OF_JOBS) for j in range(NUMBER_OF_MACHINES) for h in range(NUMBER_OF_JOBS) for g in range(NUMBER_OF_MACHINES) for k in range(NUMBER_OF_MACHINES)])

    # Assign the objective and constraint values
    out["F"] = np.array([f1, f2, f3, f4])
    out["G"] = np.array(constr_values)

algorithm = NSGA2(
pop_size=100,
sampling=BinaryRandomSampling(),
crossover=TwoPointCrossover(prob=0.9),
mutation=BitflipMutation(prob=0.1),
eliminate_duplicates=True
)

Problem initialization

problem = MyProblem()

Optimization

res = minimize(problem,
algorithm,
('n_gen', 500),
verbose=True)

Extracting the results

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))