1问题描述
某厂每日8小时的产量不低于1800件．为了进行质量控制，计划聘请两种不同水平的检验员．且每种检验员的日产量不高于1800件．一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时．检验员每错检一次，工厂要损失2元．为使总检验费用最省，该工厂应聘一级、二级检验员各几名？
2 分析问题
设需要一级和二级检验员的人数分别是x1 x2，则一天要付给检验员的费用是：
8*4*x1+8*3*x2=32x1+24x2
一天中因为检验员失误所造成的损失是：
(8*25*0.02*x1+8*15*0.05*x2)*2=8x1+12x2
约束条件：
8*25*x1+8*15*x2>=1800
8*25*x1<=1800
8*15*x2<=1800
X1>=0,x2>=0
即：
3 算法设计
li
粒子群初始化为二维数组，种群大小的行，2为列数，由所求的自变量为X1 和X2所得。
同样，粒子速度也为二维数组
粒子群算法执行过程中，加入约束条件检验函数入下：
代码：python语言
_date_ = '2019/5/27 20:16'
# -*- coding:utf-8 -*-
import numpy as np
import time
import matplotlib.pyplot as plt
import math
start = time.clock()
def getweight():
# 惯性权重
weight = 1
return weight
def getlearningrate():
# 分别是粒子的个体和社会的学习因子，也称为加速常数
# lr = (0.1,1.2)
lr = (0.49445, 1.49445)
return lr
def getmaxgen():
# 最大迭代次数
maxgen = 20
return maxgen
def getsizepop():
# 种群规模
sizepop = 20
return sizepop
def getrangepop(i):
# 粒子的位置的范围限制,x、y方向的限制相同
# rangepop = (-2*math.pi , 2*math.pi)
#rangepop = (-2,2)
if (i ==0):
rangepop = (0, 9)
else:
rangepop = (0, 15)
return rangepop
def getrangespeed():
# 粒子的速度范围限制
# rangespeed = (-0.5,0.5)
rangespeed = (-5, 5)
return rangespeed
def func(x):
# print("x的值和func")
# print(x)
if(issuitable(x)==0):
y=10000
return y
# x输入粒子位置
# y 粒子适应度值,求得最小函数值
y =40*x[0]+36*x[1]
# print(y)
return y
#是否符合约束条件
def issuitable(x):
# print("检验函数")
if(5*x[0]+3*x[1]>=45):
if((x[0]<=9.0)&(x[0]>=0.0)):
if ((x[1] <=15.0)&( x[1] >= 0.0)):
return 1
return 0
def initpopvfit(sizepop):
pop = np.zeros((sizepop,2))
v = np.zeros((sizepop,2))
fitness = np.zeros(sizepop)
for i in range(sizepop):
flag = 0
while True:
for j in range(2):
pop[i][j] = (np.random.rand() - 0.5) * getrangepop(j)[1]*2
pop[i][j] = round(pop[i][j])
pop[pop < getrangepop(j)[0]] = getrangepop(j)[0]
pop[pop > getrangepop(j)[1]] = getrangepop(j)[1]
v[i][j] = np.random.rand()-0.5
# 计算这个粒子的适应度
flag = issuitable(pop[i])
if flag == 1:
fitness[i] = func(pop[i])
break
# print("初始fitness")
# print(fitness)
return pop,v,fitness
def getinitbest(fitness,pop):
# 群体最优的粒子位置及其适应度值
gbestpop,gbestfitness = pop[fitness.argmin()].copy(),fitness.min()
#个体最优的粒子位置及其适应度值,使用copy()使得对pop的改变不影响pbestpop，pbestfitness类似
pbestpop,pbestfitness = pop.copy(),fitness.copy()
# print("全局最优")
# print(gbestpop)
# print(gbestfitness)
return gbestpop,gbestfitness,pbestpop,pbestfitness
w = getweight()
# 获取惯性权重
lr = getlearningrate()
#分别是粒子的个体和社会的学习率
maxgen = getmaxgen()
# 获得最大迭代次数
sizepop = getsizepop()
# 获取种群规模
# rangepop = getrangepop()
# 粒子的位置的范围限制,x、y方向的限制相同
rangespeed = getrangespeed()
# 粒子的速度范围限制
pop,v,fitness = initpopvfit(sizepop)
#初始化种群
gbestpop,gbestfitness,pbestpop,pbestfitness = getinitbest(fitness,pop)
# 群体最优的粒子位置及其适应度值，个体最优的粒子位置及其适应度值
result = np.zeros(maxgen)
before = np.zeros(2)
beforev = np.zeros(2)
#记录粒子的自变量的值
#返回来一个maxgen大小的一维数组；
for u in range(maxgen):
for i in range(sizepop):
flag = 0
for o in range(2):
before[o] = pop[i][o]
beforev=pop[i][o]
while True:
for j in range(2):
v[i][j] = v[i][j]+lr[0] * np.random.rand() * (pbestpop[i][j] - pop[i][j]) + lr[1] * np.random.rand() * (gbestpop[j] - pop[i][j])
v[v < rangespeed[0]] = rangespeed[0]
v[v > rangespeed[1]] = rangespeed[1]
# 粒子位置更新,需要检查粒子的约束条件是否符合
pop[i][j] = v[i][j] + pop[i][j]
pop[i][j]=round(pop[i][j])
# int是向下取整
pop[pop < getrangepop(j)[0]] = getrangepop(j)[0]
pop[pop > getrangepop(j)[1]] = getrangepop(j)[1]
flag = issuitable(pop[i])
if flag == 1:
break
else:
pop[i] = before
v[i]=beforev
#适应度更新
for j in range(sizepop):
fitness[j] = func(pop[j])
for j in range(sizepop):
if ((fitness[j] < pbestfitness[j])):
pbestfitness[j] = fitness[j]
pbestpop[j] = pop[j].copy()
if(pbestfitness.min() < gbestfitness):
gbestfitness = pbestfitness.min()
gbestpop = pop[pbestfitness.argmin()].copy()
result[u] = gbestfitness
print(result)
plt.plot(result)
plt.show()
end = time.clock()
print("time")
print((end-start)*1000.0)
由于自变量的取值空间较小，可用穷举法进行比较：
穷举法代码：
_date_ = '2019/5/28 7:43'
#利用穷举计算整数规划的值
import numpy as np
import time
start = time.clock()
num=0
max=0
m=0
n=0
b = np.zeros((160,3))
for x1 in range(10):
for x2 in range(16):
if (5 * x1 + 3 * x2>= 45):
b[num][0]=x1
b[num][1] = x2
b[num][2] = 40*x1+36*x2
num+=1
print("num")
print(num)
min=b[0][2]
m=b[0][0]
n=b[0][1]
for i in range(num):
# print("b")
# print(b)
# print("b[num][2]")
# print(b[0][2])
if(b[i][2]<min):
min=b[i][2]
m=b[i][0]
n=b[i][1]
print(min)
print(min)
print(m)
print(n)
end = time.clock()
print("time")
print((end-start)*1000.0)
#ms为单位
对比结果：
用时竟然比穷举法还长好多，最后，尝试失败
可以找到使得粒子位置更符合约束条件的方法效果应该会好很多
出现的问题：
1 速度区间问题
由用户设定用来限制粒子的速度
区间太小加上约束条件容易陷入死循环
2 约束条件问题
约束条件太复杂的话，使用随机数更新的粒子位置很难满足约束条件