Python實作粒子群（PSO）帶懲罰函數多目标優化

大佬們幫我看看，怎麼結果都一樣

import numpy as np

#--------------------------initialization_input--------------------------
Q = 15
H = 20
#--------------------------coefficients of pso--------------------------
dim = 3
wmax = 0.9
wmin = 0.4
c1 = 1.3
c2 = 1.3
itermax = 200
numpop = 20
sigma1 = 100
sigma2 = 10
M = 2
echo = 'on'
lb = [1,0.5,Q/M]
ub = [2,1,Q]
lb_array = np.tile(lb, (numpop, 1))
ub_array = np.tile(ub, (numpop, 1))
r1 = np.random.random((numpop,dim))
r2 = np.random.random((numpop,dim))
# --------------------------initialization--------------------------
t_rand = np.random.random((numpop, dim))

x = np.multiply(t_rand, (ub_array - lb_array)) + lb_array
v = np.zeros([numpop, dim], dtype=float)
vmax = 0.1 * (ub_array - lb_array)
row = np.where(x[...,0]<1.5)
x[row,0] = lb_array[row,0]
x[row,2] = Q
row = np.where(x[:,0]>=1.5)
x[row,0] = ub_array[row,0]
x[row,2] = Q/2
# --------------------------initialization calculation--------------------------
Hi = (-0.01712*x[:,2]**2 + 0.07864*x[:,2]*x[:,1] + 40.4421*x[:,1]**2).reshape(20,1)
print(Hi)
Pi = (-1.4286*0.0001*x[:,2]**3 + 0.00618*x[:,2]**2*x[:,1] + 0.04416*x[:,2]*x[:,1]**2 + 0.4402*x[:,1]**3).reshape(20,1)
# 罰函數定義
dita_max = 20
Qr_BEP = 30
Qk_BEP = (x[:,1]*Qr_BEP).reshape(20,1)
dita = (abs((x[:,2]-Qk_BEP)/Qk_BEP))
ditaH = (Hi-H).reshape(20,1)
# 揚程罰函數
[row, col] = np.where(ditaH>=0)
ditaH[row, col] = 0
[row, col] = np.where(ditaH<0)
ditaH[row, col] = ditaH[row, col]**2
# 可靠性罰函數
[row, col] = np.where(abs(dita) <= dita_max)
dita[row, col] = 0
[row, col] = np.where(abs(dita) > dita_max)
dita[row, col] = dita_max-dita[row, col]
# 目标函數求最小值
f_object1 = (-1.4286*0.0001*x[:,2].reshape(20,1)**3 +
             0.00618*x[:,2].reshape(20,1)**2*x[:,1].reshape(20,1) +
             0.04416*x[:,2].reshape(20,1)*x[:,1].reshape(20,1)**2 +
             0.4402*x[:,1].reshape(20,1)**3)*x[:,0].reshape(20,1)+\
                                                sigma1*ditaH + sigma2*dita
fbest = np.zeros([itermax,1],dtype=float)
xbest = np.zeros([itermax,3],dtype=float)
f = f_object1
fmin = np.min(f)
fmin_index = np.where(f == np.min(fmin))[0][0]
xgbest = np.multiply(np.ones([numpop, 1]), x[fmin_index, ...])
fgbest = np.multiply(np.ones([numpop, 1]), fmin)
fpbest = f
xpbest = x
j = 0
fbest[j] = fgbest[0]
xbest[j] = xgbest[0,:]

while j < itermax - 1:
    j = j + 1
    # -------update coefficients, velocity and position of pso
    w = wmax + (wmin - wmax) * j / itermax
    v = w * v + c1 * np.multiply(np.random.random((numpop, dim)), (xpbest - x)) + \
        c2 * np.multiply(np.random.random((numpop, dim)), (xgbest - x))
    x = x + v
    # -------constrain the velocity and position
    row = np.where(x[:, 0] < 1.5)
    x[row, 0] = lb_array[row, 0]
    x[row, 2] = Q
    row = np.where(x[:, 0] >= 1.5)
    x[row, 0] = ub_array[row, 0]
    x[row, 2] = Q / 2

    [row, col] = np.where(v > vmax)
    v[row, col] = vmax[row, col]
    [row, col] = np.where(v < -vmax)
    v[row, col] = -vmax[row, col]
    # -------update the fpbest and fgbest of pso
    # 罰函數定義
    Qk_BEP = (x[:, 1] * Qr_BEP).reshape(20, 1)
    dita = (abs((x[:, 2] - Qk_BEP) / Qk_BEP))
    ditaH = (Hi - H)
    # 揚程罰函數
    [row, col] = np.where(ditaH >= 0)
    ditaH[row, col] = 0
    [row, col] = np.where(ditaH < 0)
    ditaH[row, col] = ditaH[row, col] ** 2
    # 可靠性罰函數
    [row, col] = np.where(abs(dita) <= dita_max)
    dita[row, col] = 0
    [row, col] = np.where(abs(dita) > dita_max)
    dita[row, col] = dita_max - dita[row, col]

    f = f_object1
    row_min = np.where(f < fpbest)
    xpbest[row_min,:] = x[row_min,:]

    fmin = np.min(f)
    fminindex = np.where(f == np.min(fmin))[0][0]
    if fmin < fgbest[0]:
        xgbest = np.multiply(np.ones([numpop, 1]), x[fminindex, ...])
        fgbest = np.multiply(np.ones([numpop, 1]), fmin)
    fbest[j] = fgbest[0]
    xbest[j] = xgbest[0, :]


print(xbest)
           

結果

[[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]

[2. 0.74865334 7.5 ]