pso的简单实现

学了PSO（粒子群优化算法）之后，一直想自己实现一下，正好有一道例题非常合适：
已知运费和距离与货量是成正比的。某公司想要在一堆城市之间建立物流中心，物流中心和各个城市直连，并对各个城市发货。如果此公司想要节省运费成本的话，应该在何处建立此物流中心？
由于题中未能给出明确的货量大小以及城市位置，故由程序随机给出。以下是我用pso算法优化的代码：
`‘’’
我做一个迷你版的PSO算法，不知道能不能跑起来
例题背景：已知m个城市分布在二维平面中，每个城市需要一定量（不同数量）的货物，你需要建立一个仓库，
使得仓库给各个城市运送货物的运费最少，已知运费与货量和距离成正比
‘’’
import random
import math
from math import sqrt

Xmax = [20, 20] # 搜索空间
m = 10 # 城市数量
d = len(Xmax) # 维度
n = 10 # 种群数量
w = 0.7 # 权重
c1 = 2
c2 = 2
V = [] # 各粒子在各个维度的速度，初始绝对值不超过Xmax的十分之一
X = [] # 各粒子在各个空间的位置
pbest = [] # 局部最优解
gbest = [] # 全局最优解
pbestf = []
gbestf = 0.0

for i in range(n): # 初始化V，X
V.append([])
for j in range(d):
V[i].append(random.uniform(-0.1 * Xmax[j], 0.1 * Xmax[j]))
X.append([])
for j in range(d):
X[i].append(random.uniform(0, Xmax[j]))

def randstart(m, Xmax): # 初始随机条件
np = []
for i in range(m):
np.append(random.uniform(0, 10)) # 随机范围自定，这里选取最大为10
Xp = []
for i in range(m):
Xp.append([])
for j in range(d):
Xp[i].append(random.uniform(0, Xmax[j]))
return np, Xp # 这里需要两个变量来接收函数返回值，一个是各个城市的货量，一个是各个城市的位置

np, Xp = randstart(m, Xmax) # 初始化np，Xp

def f(x): # 适应值判断,此时f（x）越小越优
sum = 0
for i in range(m):
L = 0
for j in range(d):
l = math.pow(Xp[i][j] - x[j], 2)
L += l
sum += np[i] * sqrt(l)
return sum / 300

def ran(): # 返回（0,1）随机数
return random.uniform(0, 1)

for i in range(n): # pbest，gbest初始化
pbest.append([])
for j in range(d):
pbest[i].append(0)
pbest[i][j] = X[i][j]
if i == 0:
gbest.append(0)
gbest[j] = X[i][j]
pbestf.append(0.0)
pbestf[i] = f(X[i])
if i == 0: # gbestf初始化
gbestf = f(X[i])
else:
if pbestf[i] < gbestf:
gbestf = pbestf[i]
for j in range(d):
gbest[j] = pbest[i][j]
#print(pbest)

for t in range(10): # 开始进行粒子群搜索，t为搜索次数
for i in range(n):
for j in range(d):
V[i][j] = w * V[i][j] + c1 * ran() * (pbest[i][j] - X[i][j])
+ c2 * ran() * (gbest[j] - X[i][j]) # 速度更新公式
X[i][j] = X[i][j] + V[i][j] # 位置更新公式，主要这边位置有越过边界的可能，但是例题不需要考虑此种情况
for i in range(n):
fx = f(X[i]) # 防止重复运算f(X[i])
if fx < pbestf[i]: # 更新pbest
pbestf[i] = fx
for j in range(d):
pbest[i][j] = X[i][j]
if pbestf[i] < gbestf: # 更新gbest
gbestf = pbestf[i]
for j in range(d):
gbest[j] = pbest[i][j]
print(gbestf)
print(gbest)`

据说比较简单的方法是，选择一个轻薄木板，然后对应城市位置上穿孔，各个空上按照比例加重物，最后看看此木板在哪个位置上平衡。啊，但是作为一只预备程序猿，代码什么的，毛毛雨啦~