目录
１配电网故障定位原理
2 二进制粒子群算法的简单举例------背包问题
2.1 算例
2.2 Python 实现
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm#进度条设置
import matplotlib.pyplot as plt
import matplotlib; matplotlib.use('TkAgg')
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
class BackPack(object):
def __init__(self):
self.Weight = [95,75,23,73,50,22,6,57,89,98] # 物品体积
self.Value = [89, 59, 19, 43, 100, 72, 44, 16, 7, 64] # 物品价值
self.N = 100 # 群体粒子个数
self.D = len(self.Weight) # 粒子维数
self.T = 100 # 最大迭代次数
self.c1 = 1.5 # 学习因子1
self.c2 = 1.5 # 学习因子2
self.w = 1 # 惯性因子，一般取1
self.V_max = 10 # 速度最大值
self.V_min = -10 # 速度最小值
self.Weight_max = 300 # 背包容量
self.afa = 10 # 惩罚系数
self.G=100 #迭代次数
#=====初始化种群=========
def init_x(self):
"""
:return: 随机生成的种群(二维list)
"""
X = np.random.choice([0, 1], size=(self.N, self.D))
return X
#======初始化速度==========
def init_v(self):
"""
:return: 速度
"""
V = np.random.random(size=(self.N, self.D)) ##10维
return V
#=====适应度值========
def fitness_func(self,X):
"""
:param X: 种群
:return:
"""
"""计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 10 """
size = np.shape(X)[0] #种群个数
result = np.zeros((size, 1))
for i in range(size): # 遍历每一个粒子
fit = sum(X[i] * self.Value) # 物品价值
TotalSize = sum(X[i] * self.Weight) # 物品重量
if TotalSize <= self.Weight_max: # 如果小于限制重量
fit = fit
result[i] = fit
else:
#fit = fit - self.afa * (TotalSize - self.Weight_max) # 对fit进行惩罚
fit=0
result[i] = fit
return result # 我们要求result越大
#=====速度更新公式========
def velocity_update(self,V, X, pbest, gbest):
"""
根据速度更新公式更新每个粒子的速度
种群size=100
:param V: 粒子当前的速度矩阵，size*10 的矩阵
:param X: 粒子当前的位置矩阵，size*10 的矩阵
:param pbest: 每个粒子历史最优位置，size*10 的矩阵
:param gbest: 种群历史最优位置，1*10 的矩阵
"""
r1 = np.random.random((self.N, 1)) #(粒子个数,1)
r2 = np.random.random((self.N, 1))
V = self.w * V + self.c1 * r1 * (pbest - X) + self.c2 * r2 * (gbest - X) # 直接对照公式写就好了
#===防止越界处理===
V[V < self.V_min] = np.random.random() * (self.V_max - self.V_min) + self.V_min
V[V > self.V_max] = np.random.random() * (self.V_max - self.V_min) +self.V_min
return V
#==========位置更新公式==========
def position_update(self,X, V):
"""
根据公式更新粒子的位置
:param X: 粒子当前的位置矩阵，维度是 size*10
:param V: 粒子当前的速度举着，维度是 size*10
return 更新后的种群
"""
for i in range(self.N): # 遍历每一个粒子
# 修改速度为sigmoid形式
V[i, :] = 1. / (1 + np.exp(-np.array(V[i, :])))
for j in range(self.D): # 遍历粒子中的每一个元素
rand = np.random.random() # 生成 0-1之间的随机数
if V[i, j] > rand:
X[i, j] = 1
else:
X[i, j] = 0
#对当前个体进行限制，放在超重
while np.sum(X[i]*self.Weight)>self.Weight_max:#如果当前粒子超重
X[i]=np.random.choice([0, 1], size=(1, self.D))
return X
def update_pbest(self,X, fitness, pbest, pbestfitness, m):
"""
更新个体最优
:param X: 当前种群
:param fitness: 当前每个粒子的适应度
:param pbest: 更新前的个体最优解
:param pbestfitness: 更新前的个体最优适应度
:param m: 粒子数量
:return: 更新后的个体最优解、个体最优适应度
"""
for i in range(m):
if fitness[i] > pbestfitness[i]:
pbest[i] = X[i]
pbestfitness[i] = fitness[i]
return pbest, pbestfitness
def update_gbest(self,pbest, pbestfitness, gbest, gbestfitness, m):
"""
更新全局最优解
:param pbest: 粒子群
:param pbestfitness: 个体适应度(个体最优)
:param gbest: 全局最优解
:param gbestfitness: 全局最优解
:param m: 粒子数量
:return: gbest全局最优，g对应的全局最优解
"""
for i in range(m):
if pbestfitness[i] > gbestfitness:
gbest = pbest[i]
gbestfitness = pbestfitness[i]
return gbest, gbestfitness
def main(self):
fitneess_value_list = [] # 记录每次迭代过程中的种群适应度值变化
x = self.init_x() # 初始化x
v = self.init_v() # 初始化v
# 计算种群各个粒子的初始适应度值
p_fitness = self.fitness_func(x)
# 计算种群的初始最优适应度值
g_fitness = p_fitness.max()
# 讲添加到记录中
fitneess_value_list.append(g_fitness)
# 初始的个体最优位置和种群最优位置
pbest = x
gbest = x[p_fitness.argmax()] #
#===接下来就是不断迭代了========
for i in tqdm(range(self.G)):
pbest=pbest.copy()#必须加，不然会被篡改值，造成结果错
p_fitness= p_fitness.copy()#必须加，不然会被篡改值，造成结果错
gbest=gbest.copy()#必须加，不然会被篡改值，造成结果错
g_fitness=g_fitness.copy()#必须加，不然会被篡改值，造成结果错
v = self.velocity_update(v, x, pbest, gbest) # 更新速度
x = self.position_update(x, v) # 更新位置
p_fitness2 = self.fitness_func(x) # 计算子代各个粒子的适应度
# 更新每个粒子的历史最优位置
pbest, p_fitness=self.update_pbest(x, p_fitness2, pbest, p_fitness, self.N)
#更新群体的最优位置
gbest, g_fitness=self.update_gbest(pbest, p_fitness, gbest, g_fitness, self.N)
# 记录最优迭代结果
fitneess_value_list.append(g_fitness)
print("最优适应度是：%.5f" % fitneess_value_list[-1])
print("最优解是", gbest)
print('最优解对应的体积',np.sum(gbest*self.Weight))
print('最优解对应的价值',np.sum(gbest*self.Value))
plt.plot(fitneess_value_list,label='迭代曲线')
plt.xlabel('迭代次数')
plt.ylabel('适应度')
plt.legend()
plt.show()
if __name__=='__main__':
pso = BackPack()
pso.main()
2.3 Matlab 实现
%% 离散粒子群算法解决0-1背包问题
%% 初始化
clear all; %清除所有变量
close all; %清图
clc; %清屏
N=100; %群体粒子个数
D=10; %粒子维数
c1=1.5; %学习因子1
c2=1.5; %学习因子2
Wmax=0.8; %惯性权重最大值
Wmin=0.4; %惯性权重最小值
Vmax=10; %速度最大值
Vmin=-10; %速度最小值
V = 300; %背包容量
C = [95,75,23,73,50,22,6,57,89,98]; %物品体积
W = [89,59,19,43,100,72,44,16,7,64]; %物品价值
afa = 2; %惩罚函数系数
%% 初始化种群个体（限定位置和速度）
x=randi(N,D); %随机获得二进制编码的初始种群，randint在后期版本中会被randi代替
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%% 初始化个体最优位置和最优值
p=x;
pbest=ones(D,1);
for i=1:D
pbest(i)= func(x(i,:),C,W,V,afa);
end
%% 初始化全局最优位置和最优值
g=ones(1,D);
gbest=-10000;
for i=1:D
if(pbest(i)>gbest)
g=p(i,:);
gbest=pbest(i);
end
end
gb=ones(1,D);
%% 按照公式依次迭代直到满足精度或者迭代次数
while 1>0 %无限迭代
for j=1:D
%% 更新个体最优位置和最优值
if (func(x(j,:),C,W,V,afa)>pbest(j))
p(j,:)=x(j,:);
pbest(j)=func(x(j,:),C,W,V,afa);
end
%% 更新全局最优位置和最优值
if(pbest(j)>gbest)
g=p(j,:);
gbest=pbest(j);
end
%% 计算动态惯性权重值
w=Wmax-(Wmax-Wmin)*i/D;
%% 跟新位置和速度值
v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
+c2*rand*(g-x(j,:));
%% 边界条件处理
for ii=1:D
if (v(j,ii)>Vmax) | (v(j,ii)< Vmin)
v(j,ii)=rand * (Vmax-Vmin)+Vmin;
end
end
vx(j,:)=1./(1+exp(-v(j,:)));
for jj=1:D
if vx(j,jj)>rand
x(j,jj)=1;
else
x(j,jj)=0;
end
end
end
%% 记录历代全局最优值
gb(i)=gbest;
if g.*C<=300 %如果满足条件 退出
break
end
g; %最优个体
figure
plot(gb)
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线')
end
%% 适应度函数
function result = func(f,C,W,V,afa)
fit = sum(f.*W);
TotalSize = sum(f.*C);
if TotalSize <= V
fit = fit;
else
fit = fit - afa * (TotalSize - V);
end
result = fit;3 二进制粒子群算法的配电网故障定位(Matlab实现)
clc;
clear all;
%% 初始化参数
global Ij; %全局变量
global Vmax;%允许误差
MaxNum=120;%粒子最大迭代次数
D=12;%粒子的维数代表配电网的馈线区段总数
Particlesize=30;%粒子群规模
c1=2.1;%每个粒子的个体学习因子，也称为加速常数
c2=2.1;%每个粒子的社会学习因子，也称为加速常数
w=1;%惯性因子
Vmax=4;%粒子的最大飞翔速度
S_b=randi([0,1],Particlesize,D);%粒子所在的位置,randint在后期版本中会被randi代替
V=Vmax*rand(Particlesize,D);%粒子的飞翔速度
Ij=[1 1 1 1 0 1 1 1 0 0 1 1];%FTU 反馈的值这里就是输入，请修改这里，看结果
fitness=F(S_b');%这里注意下，要转置
%% inline定义的适应度函数会使程序运行速度大大降低
for i=1:Particlesize
f(i)=F(S_b(i,:));%S_b配电网中各设备状态,为1表明设备故障,取0则设备正常。
end
%% 个体最优与总群最优
personalbest_x=S_b;%这里是种群里的每个粒子都认为自己是最好的
personalbest_faval=f;
[globalbest_favali]=min(personalbest_faval);%全局最优
globalbest_x=personalbest_x(i,:);
k=1;
Start=0;
%% 这里就是把整个种群中的所有粒子都检查了一遍，保证当前最优
while k<=MaxNum
Remember=globalbest_favali;
for i=1:Particlesize
f(i)=F(S_b(i,:));%计算适应度函数值
if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置
personalbest_faval(i)=f(i);
personalbest_x(i,:)=S_b(i,:);
end
end
%% 更新种群最优
[globalbest_favali]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
%% 更新粒子群里每个个体的最新位置
for i=1:Particlesize
V(i,:)=w*V(i,:)+c1*rand*(personalbest_x(i,:)-S_b(i,:))...
+c2*rand*(globalbest_x-S_b(i,:));%更新了速度
for j=1:length(Ij)
if rand(1)<Sigmoid(V(i,j))%这里体现了速度决定位置
S_b(i,j)=1;
else
S_b(i,j)=0;
end
end
end
%% 粒子的位置代表配电网中馈线区段的状态globalbest_x
if globalbest_favali-Remember==0
Start=Start+1;
if Start==20
disp(globalbest_x);
break;
end
end
k=k+1;
end
function fitness=F(SB)
global Ij;
%% 期望函数Ij(SB):IE
IE=zeros(1,length(Ij));
IE(1)=SB(1)|SB(2)|SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);
IE(2)=SB(2)|SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);
IE(3)=SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);
IE(4)=SB(4)|SB(5);
IE(5)=SB(5);
IE(6)=SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);
IE(7)=SB(7)|SB(8)|SB(9)|SB(10);
IE(8)=SB(8)|SB(9)|SB(10);
IE(9)=SB(9);
IE(10)=SB(10);
IE(11)=SB(11)|SB(12);
IE(12)=SB(12);
%% 构造配电网故障定位评价函数Fij:fitness
fitness=sum(abs(Ij-IE))+0.5*sum(abs(SB));
end
function AN=Sigmoid(V)
global Vmax;
if V>Vmax
AN=0.98;
elseif V>-Vmax
AN=1/(1+exp(-V));
else
AN=-0.98;
end
end 
