'''
Created on 2019年5月24日
@author: 紫薇星君
'''
import pandas as pd
from sklearn import preprocessing,model_selection
import numpy as np
from sklearn import svm
from matplotlib import pyplot as plt
import pickle
from mpl_toolkits.mplot3d import Axes3D
class F():
def __init__(self,x):
self.x=x
self.mi=x.min()
self.ma=x.max()
def f(self,a):#归一化
return (a-self.mi)/(self.ma-self.mi)
def f_b(self,b):#反归一化
return b*(self.ma-self.mi)+self.mi
def readpkl(road):#存储的数据读取函数
file_pkl=open(road,'rb')
readed=pickle.load(file_pkl)
file_pkl.close()
return readed
def savepkl(road,data_pkl):#pkl文件存储函数
file_pkl=open(road,'wb')
pickle.dump(data_pkl,file_pkl)
file_pkl.close()
return
class PSO():
'''参数设定'''
def __init__(self):
self.w = self.getweight()
self.lr = self.getlearningrate()
self.maxgen = self.getmaxgen()
self.sizepop = self.getsizepop()
self.rangepop = self.getrangepop()
self.rangespeed = self.getrangespeed()
return
def getweight(self):
# 惯性权重
weight = 0.5
return weight
def getlearningrate(self):
# 分别是粒子的个体和社会的学习因子，也称为加速常数
lr = (0.49445,1.49445)
return lr
def getmaxgen(self):
# 最大迭代次数
maxgen = 500
return maxgen
def getsizepop(self):
# 种群规模
sizepop = 200
return sizepop
def getrangepop(self):
# 粒子的位置的范围限制,x、y方向的限制相同
rangepop = (-1,1)
return rangepop
def getrangespeed(self):
# 粒子的速度范围限制
rangespeed = (-0.001,0.001)
return rangespeed
'''end'''
'''种群初始化'''
def func(self,x):
# x输入粒子位置
# y 粒子适应度值
x1=[x[0],x[1],t_set]
y = fpn.f_b(clf_svr.predict([x1]))
return y
'''然后初始化粒子和计算适应度值'''
def initpopvfit(self,sizepop):
pop = np.zeros((sizepop,2))#所有鸟的位置矩阵，几个输入就有几维
v = np.zeros((sizepop,2))#所有鸟的速度矩阵，几个输入就有几个分速度，就有几维
fitness = np.zeros(sizepop)#所有鸟的适应度列表
for i in range(sizepop):
pop[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]#鸟位置初始化
v[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]#鸟速度初始化
fitness[i] = self.func(pop[i])#适应度函数初始化
return pop,v,fitness
'''end'''
'''寻找初始化后的鸟最优位置与适应度值'''
def getinitbest(self,fitness,pop):
# 群体最优的粒子位置及其适应度值
gbestpop,gbestfitness = pop[fitness.argmax()].copy(),fitness.max()#a.argmax()指数组a最大值的位置;a.max()指a的最大值
#个体最优的粒子位置及其适应度值,使用copy()使得对pop的改变不影响pbestpop，pbestfitness类似
pbestpop,pbestfitness = pop.copy(),fitness.copy()
return gbestpop,gbestfitness,pbestpop,pbestfitness
'''end'''
'''迭代寻优'''
def run(self):
pop, v, fitness = self.initpopvfit(self.sizepop)
gbestpop, gbestfitness, pbestpop, pbestfitness = self.getinitbest(fitness, pop)
result = np.zeros(self.maxgen)#建立一个装每次迭代运算结果的数组
for i in range(self.maxgen):#迭代计算
# 速度更新
for j in range(self.sizepop):#遍历所有鸟
v[j] += self.lr[0] * np.random.rand() * (pbestpop[j] - pop[j]) + self.lr[1] * np.random.rand() * (
gbestpop - pop[j])#Ir是非负加速因子，np.random.rand()是(0,1)上的随机数，防止鸟满目搜索，得到新速度矩阵
v[v < self.rangespeed[0]] = self.rangespeed[0]#rangespeed[0]是速度下限,将v中所有小于下限的值全改为下限
v[v > self.rangespeed[1]] = self.rangespeed[1]#rangespeed[1]是速度上限，将v中所有大于上限的值全改为上限
# 粒子位置更新
for j in range(self.sizepop):#遍历所有鸟
pop[j] += 0.5 * v[j]#得到新的位置列表
pop[pop < self.rangepop[0]] = self.rangepop[0]#鸟位置的下限，即搜索的范围最小值点
pop[pop > self.rangepop[1]] = self.rangepop[1]#鸟位置的上限，即搜索的范围最大值点
# 适应度更新
for j in range(self.sizepop):#遍历所有鸟
fitness[j] = self.func(pop[j])#新的适应度列表
for j in range(self.sizepop):#遍历所有鸟
if fitness[j] > pbestfitness[j]:
pbestfitness[j] = fitness[j]#所有鸟的各自经历的最大适应度值
pbestpop[j] = pop[j].copy()#上面最大适应度值时对应的位置
if pbestfitness.max() > gbestfitness:
gbestfitness = pbestfitness.max()#所有鸟经历过的最大适应度值
gbestpop = pop[pbestfitness.argmax()].copy()#该适应度对应的位置
result[i] = gbestfitness#此次迭代寻优的结果
return result,gbestfitness,gbestpop
'''end'''
if __name__=='__main__':
'''归一输入：[[t1,li1,pfd1],[t2,li2,pfd2]...] 归一输出：[[pn1],[pn2]...]'''
road='.\datafile\jiaohuan2.csv'
file=pd.read_csv(road)
t=np.array(file['T'])
lq=np.array(file['LQ'])
pn=np.array(file['Pn'])
pfd=np.array(file['PFD'])
num=len(pn)
'''*************表格处理结束**************'''
'''归一化'''
ft=F(t)
flq=F(lq)
fpn=F(pn)
fpfd=F(pfd)
t_set1=32
t_set=ft.f(t_set1)
clf_svr=readpkl('.\savefile\clf_pnflour.pkl')
pso = PSO()
progress,result,point = pso.run()
plt.figure()
plt.plot(progress)
plt.xlabel(u'迭代次数',fontproperties='SimHei')
plt.ylabel(u'净光合速率/(μmol/(m$^2$·s))',fontproperties='SimHei')
plt.title('温度'+':'+'T'+'='+str(t_set1)+'℃',fontproperties='SimHei')
plt.show()
print(result)
point_1=(t_set,point[0],point[1])
t_set=ft.f_b(t_set)
point_2=[fpfd.f_b(point[0]),flq.f_b(point[1]),t_set]
print(point_2)
'''设置当前数据'''
pfd_set=np.array([i*3 for i in range(600)])
lq_set=np.array([i*0.02 for i in range(50)])
t_1=ft.f(t_set)
pfd_1=fpfd.f(pfd_set)
lq_1=flq.f(lq_set)
'''计算'''
pn=np.empty((len(lq_1),len(pfd_1)))
for i in range(len(lq_1)):
for j in range(len(pfd_1)):
pn[i][j]=fpn.f_b(clf_svr.predict([[pfd_1[j],lq_1[i],t_1]]))
'''画图'''
Y,X=np.meshgrid(pfd_set,lq_set)
Z=pn
fig=plt.figure()
ax=Axes3D(fig)
ax.set_xlabel('Li')
ax.set_ylabel('PFD')
ax.set_zlabel('Pn')
ax.plot_surface(X,Y,Z)
plt.title('T='+str(t_set))
ax.scatter([point_2[1]],[point_2[0]],[result],marker='*',c='r')
plt.show()
'''结束'''
评论0