adding ensure constraint for Particle Swarm Optimization

Question

I tried to add some ensure constraints to make sure that my PSO result satisfied the constraint. However when i run the code, it does not really work. My input data: https://docs.google.com/spreadsheets/d/16F4wvrUXYRMl5KbUiFvSYoLDwGg9qlxs/edit?usp=sharing&ouid=108297438086372368891&rtpof=true&sd=true

Here is my python code: import numpy as np import pandas as pd from pulp import * import random as rd

rd.seed(15)
W = 0.5
c1 =0.5
c2 =0.5
n_iterations = 5
n_particles = 5
target_error = 1e-6

df = pd.read_excel(r'D:\final\DATA.xlsx', sheet_name='Set1', header=None)
v = np.array(df.iloc[1:4, 2:6].values.tolist(), dtype=float)
f = np.array(df.iloc[6:9, 2:6].values.tolist(), dtype=int)
u1 = np.array(df.iloc[12:16, 2:6].values.tolist(), dtype=float)
u2 = np.array(df.iloc[12:16, 9:13].values.tolist(), dtype=float)
u3 = np.array(df.iloc[12:16, 15:19].values.tolist(), dtype=float)
C = np.array(df.iloc[20:23, 2:6].values.tolist(), dtype=float)
MOQ = np.array(df.iloc[26:29, 2:6].values.tolist(), dtype=float)
D = np.array(df.iloc[41, 1:4].values.tolist(), dtype=int)
strat = np.array(df.iloc[32:35, 2:6].values.tolist(), dtype=int)
Nmax = np.array(df.iloc[37, 1:4].values.tolist(), dtype=int)
Nstr = np.array(df.iloc[37, 5:8].values.tolist(), dtype=int)
P = range(0,2,1)
I = range(0,3,1)
print(I, P)
u = {}
    for i in I:
    for j in I:
        for p in P:
            if p == 0:
                u[i, j, p] = u1[i, j]
            elif p == 1:
                u[i, j, p] = u2[i, j]
            else:
                u[i, j, p] = u3[i, j]


class Particle():
    def __init__(self):
        y = np.array([[1, 0, 1, 0],
                  [0, 1, 1, 0],
                  [1, 0, 0, 1]], dtype=int)
        x = np.array([[0.5, 0, 0.5, 0],
                  [0, 0.6, 0.4, 0],
                  [0.4, 0, 0, 0.6]], dtype=float)
        a = round(np.random.uniform(-0.2, 0.2),2)
        b = round(np.random.uniform(-0.2, 0.2),2)
        c = round(np.random.uniform(-0.2, 0.2),2)

        offset = np.array([[a, 0, -a, 0], [0, b, -b, 0], [c, 0, 0, -c]])
        x += offset

        self.position = np.array([[x], [y]])  # Matrix form
        self.pBest_position = self.position
        self.pBest_value = float('inf')
        self.velocity = np.zeros(self.position.shape)

    def update(self):
        self.position = self.position + self.velocity

class Space():
    def __init__(self, target, target_error, n_particles):
        self.target = target
        self.target_error = target_error
        self.n_particles = n_particles
        self.particles = []
        self.gBest_value = float('inf')
        self.gBest_position = np.random.rand(2, 1) * 50  # Matrix form

    def fitness(self, particle):
        x = particle.position[0][0]  
        y = particle.position[1][0]  
        # Calculation for Z
        Z = 0  # Initialize Z as 0
        for i in I:
            for p in P:
                Z += (D[p] * v[i][p]) * x[i][p] + f[i][p] * y[i][p]

        # Calculation for K
        K = 0  # Initialize K as 0
        for i in I:
            for j in I:
                for p in P:
                    K += lpSum(u[i, j, p] * x[i][p] * x[j][p])
        F = Z + K
        return F

    def set_pBest(self):
        for particle in self.particles:
            fitness_value = value(self.fitness(particle))  
            if particle.pBest_value > fitness_value:
                particle.pBest_value = fitness_value
                particle.pBest_position = particle.position.copy()

    def set_gBest(self):
        for particle in self.particles:
            if self.gBest_value > particle.pBest_value:
                self.gBest_value = particle.pBest_value
                self.gBest_position = particle.pBest_position.copy()

    def update_particles(self):
        for particle in self.particles:
            r1 = np.random.random(particle.position.shape)
            r2 = np.random.random(particle.position.shape)
            particle.velocity = (W * particle.velocity +
                                 c1 * r1 * (particle.pBest_position - particle.position) +
                                 c2 * r2 * (self.gBest_position - particle.position))
            particle.update()
# Ensure constraint 1: sum x[i][p] = 1 
            for i in I:
                for p in P:
                    sum_p = np.sum(particle.position[0][0][i])  # Calculate the sum for the current i     index
                    diff = 1 - sum_p  # Calculate the difference from 1
                    particle.position[0][0][i][p] += diff / len(P)  # Add the difference equally to each value for the current i index
# Ensure constraint 6 & 7:
            for i in I:
                for p in P:
                    if MOQ[i][p] * particle.position[1][0][i][p] > particle.position[0][0][i][p]:
                        particle.position[0][0][i][p] = MOQ[i][p] * particle.position[1][0][i][p]
                    elif particle.position[0][0][i][p] > C[i][p] * particle.position[1][0][i][p]:
                        particle.position[0][0][i][p] = C[i][p] * particle.position[1][0][i][p]
                
space = Space(target=0, target_error=target_error, n_particles=n_particles)

# Initialize particles
for _ in range(n_particles):
    particle = Particle()
    space.particles.append(particle)

iteration = 0
while iteration < n_iterations:
    space.set_pBest()
    space.set_gBest()

    if abs(space.gBest_value - space.target) <= space.target_error:
        break

    space.update_particles()

    # Results for each iteration
    print("Iteration:", iteration + 1)
    print("Best F:", space.gBest_value)
    print("Best x:")
    print(space.gBest_position[0][0])
    print("Best y:")
    print(space.gBest_position[1][0])
    print("------------------------")

    iteration += 1

# Final results
print("Final Best F:", space.gBest_value)
print("Final Best x:")
print(space.gBest_position[0][0])
print("Final Best y:")
print(space.gBest_position[1][0])


# Constraint: forall(p in P) sum(i in I) x[i][p] == 1
x = space.gBest_position[0][0]  # Extract x value from matrix
y = space.gBest_position[1][0]  # Extract y value from matrix

for i in I:
    if sum(x[i]) != 1:
        print("Constraint 1 not satisfied: forall(p in P) sum(i in I) x[i][p] == 1")

# Constraint: forall(p in P) (sum(i in I) (D[1][p]*v[i][p]*x[i][p] + f[i][p]*y[i][p]) <= B)
for p in P:
    if sum(D[p] * v[i][p] * x[i][p] + f[i][p] * y[i][p] for i in I) > 90000:
        print("Constraint 2 not satisfied: forall(p in P) (sum(i in I) (D[1][p]*v[i][p]*x[i][p] + f[i][p]*y[i][p]) <= 90000)")

# Constraint: forall(p in P) 2 <= sum(i in I) y[i][p]
for i in I:
    if sum(y[i]) < 2:
        print("Constraint 3 not satisfied: forall(p in P) 2 <= sum(i in I) y[i][p]")

# Constraint: forall(p in P) sum(i in I) y[i][p] <= Nmax[p]
for p in P:
    if sum(y[i][p] for i in I) > Nmax[p]:
        print("Constraint 4 not satisfied: forall(p in P) sum(i in I) y[i][p] <= Nmax[p]")

# Constraint: forall(p in P) Nstr[p] <= sum(i in I) y[i][p]*strat[i][p]
for p in P:
    if Nstr[p] > sum(y[i][p] * strat[i][p] for i in I):
        print("Constraint 5 not satisfied: forall(p in P) Nstr[p] <= sum(i in I) y[i][p]*strat[i][p]")

# Constraint: forall(i in I, p in P) x[i][p] <= C[i][p]*y[i][p]
for i in I:
    for p in P:
        if x[i][p] > C[i][p] * y[i][p]:
            print("Constraint 6 not satisfied: forall(i in I, p in P) x[i][p] <= C[i][p]*y[i][p]")

# Constraint: forall(i in I, p in P) MOQ[i][p]*y[i][p] <= x[i][p]
for i in I:
    for p in P:
        if MOQ[i][p] * y[i][p] > x[i][p]:
            print("Constraint 7 not satisfied: forall(i in I, p in P) MOQ[i][p]*y[i][p] <= x[i][p]")

# Constraint: forall(i in I, p in P) 0 <= x[i][p] <= 1, 0 <= y[i][p] <= 1, 0 <= C[i][p] <= 1, 0 <= MOQ[i][p] <= 1, 0 <= strat[i][p] <= 1
for i in I:
    for p in P:
        if not (0 <= x[i][p] <= 1 and 0 <= y[i][p] <= 1 and 0 <= C[i][p] <= 1 and 0 <= MOQ[i][p] <= 1 and 0 <= strat[i][p] <= 1):
            print("Constraint 8 not satisfied: forall(i in I, p in P) 0 <= x[i][p] <= 1, 0 <= y[i][p] <= 1, 0 <= C[i][p] <= 1, 0 <= MOQ[i][p] <= 1, 0 <= strat[i][p] <= 1")