def doy2date(year,doy): month_leapyear=[31,29,31,30,31,30,31,31,30,31,30,31] month_notleap= [31,28,31,30,31,30,31,31,30,31,30,31] if year%4==0 and (year%100!=0 or year%400==0): for i in range(0,12): if doy>month_leapyear[i]: doy-=month_leapyear[i] continue if doy<=month_leapyear[i]: month=i+1 day=doy break else: for i in range(0,12): if doy>month_notleap[i]: doy-=month_notleap[i] continue if doy<=month_notleap[i]: month=i+1 day=doy break return month,day def form_data(id): day = int((id)/48)+1 year = 2019 month,days = doy2date(year,day) hour = int((id%48)/2) mins =int((id%48)%2)*30 month_leapyear=[31,29,31,30,31,30,31,31,30,31,30,31] if days>month_leapyear[month-1]: month+=1 day=1 return str(year)+'-'+str(month)+'-'+str(days)+' '+str(hour)+':'+str(mins) def read_y(path = 'y(1).csv'): df = pd.read_csv(path) df.head() df = df.set_index('Unnamed: 0') data_array = df.values.reshape(-1) data = data_array df = pd.DataFrame(data) df = df.reset_index() df.columns = ['id','value'] df['date'] = df['id'] .apply(lambda x: form_data(x)) df['date'] = pd.to_datetime(df['date']) return df def imf_n(data, num, bc="natural"): imfs = emd(data, bc=bc) if len(imfs) >= num: return np.array(imfs[num - 1]) else: return None def ceemdan(data, I=500, sd=0.1, max_extr=2, bc="natural",Imfs_max = 5): epsilon = 0.08 size = data.shape[0] imfs = [] reside = data.copy() imfs1 = [] for i in range(int(I)): imf1 = imf_n(reside + epsilon * np.random.normal(0, 1, size), 1, bc=bc) imfs1.append(imf1) imf1 = sum(imfs1) / I imfs.append(imf1) print("imf 0 is ready") num = 1 last_imf = False for iteration in range(Imfs_max): iteration += 1 imf_prev = imfs[num - 1] reside -= imf_prev imfs_num = [] for i in range(int(I)): imf_num = imf_n(reside + epsilon * imf_n(np.random.normal(0, 1, size), num, bc=bc), 1, bc=bc) if imf_num is None: last_imf = True break imfs_num.append(imf_num) if not last_imf: imf_num = sum(imfs_num) / I # finish criterion by number of extremas if (len(signal.argrelmax(imf_num, order=1)[0]) + len(signal.argrelmin(imf_num, order=1)[0])) <= max_extr: last_imf = True # finish criterion by difference in amps sd_ = np.sum((imf_num - imf_prev) ** 2) / np.sum(imf_prev ** 2) if sd_ < sd: last_imf = True if last_imf: imfs.append(np.array(reside)) print("last imf is ready") break imfs.append(imf_num) print("imf", num, "is ready") num += 1 return imfs def waveletdec(s, coef_type='d', wname='sym7', level=6, mode='symmetric'): import pywt N = len(s) w = pywt.Wavelet(wname) a = s ca = [] cd = [] for i in range(level): (a, d) = pywt.dwt(a, w, mode) # 将a作为输入进行dwt分解 ca.append(a) cd.append(d) rec_a = [] rec_d = [] for i, coeff in enumerate(ca): coeff_list = [coeff, None] + [None] * i rec_a.append(pywt.waverec(coeff_list, w)[0:N]) # 进行重构 for i, coeff in enumerate(cd): coeff_list = [None, coeff] + [None] * i rec_d.append(pywt.waverec(coeff_list, w)[0:N]) # 进行重构 if coef_type == 'd': return rec_d return rec_a def weight_variables(shape): w = tf.Variable(tf.random_normal(shape=shape, mean=0.0, stddev=1.0)) return w def bias_variables(shape): b = tf.Variable(tf.constant(0.001, shape=shape)) return b def single_RNN(x, hidden_nodes0, hidden_nodes, input_features, output_class): x_reshape = tf.reshape(x, [-1, 1, input_features]) with tf.variable_scope("LSTM"): rnn_cell = tf.nn.rnn_cell.MultiRNNCell( [tf.nn.rnn_cell.LSTMCell(hidden_nodes0), tf.nn.rnn_cell.LSTMCell(hidden_nodes0)]) outputs, _ = tf.nn.dynamic_rnn(rnn_cell, x_reshape, dtype=tf.float32) x_at_reshape = tf.reshape(outputs, [-1, hidden_nodes0]) # 全连接分类层 with tf.variable_scope("OUTPUT"): w_fc2 = weight_variables([hidden_nodes0, output_class]) b_fc2 = bias_variables([output_class]) y_predict = tf.matmul(x_at_reshape, w_fc2) + b_fc2 return tf.reshape(y_predict, [-1, output_class]) # 双向LSTM def RNN(x, hidden_nodes0, hidden_nodes, input_features, output_class): x_reshape = tf.reshape(x, [-1, 1, input_features]) # 双向lstm层 with tf.variable_scope("BILSTM"): rnn_cellforword = tf.nn.rnn_cell.MultiRNNCell( [tf.nn.rnn_cell.LSTMCell(hidden_nodes0), tf.nn.rnn_cell.LSTMCell(hidden_nodes0)]) rnn_cellbackword = tf.nn.rnn_cell.MultiRNNCell( [tf.nn.rnn_cell.LSTMCell(hidden_nodes), tf.nn.rnn_cell.LSTMCell(hidden_nodes)]) outputs, _ = tf.nn.bidirectional_dynamic_rnn(rnn_cellforword, rnn_cellbackword, x_reshape, dtype=tf.float32) nodes = hidden_nodes0 + hidden_nodes outputs = tf.concat(outputs, 2) x_at_reshape = tf.reshape(outputs, [-1, nodes]) # 全连接分类层 with tf.variable_scope("OUTPUT"): w_fc2 = weight_variables([nodes, output_class]) b_fc2 = bias_variables([output_class]) y_predict = tf.matmul(x_at_reshape, w_fc2) + b_fc2 return tf.reshape(y_predict, [-1, output_class]) def fun(pop,P,T,Pt,Tt): tf.reset_default_graph() tf.set_random_seed(0) num_epochs = int(pop[0])#迭代次数 alpha = pop[1]# 学习率 hidden_nodes0 = int(pop[2])#第一隐含层神经元 hidden_nodes = int(pop[3])#第二隐含层神经元 input_features = P.shape[1] output_class = T.shape[1] batch_size = 128# batchsize # placeholder X = tf.placeholder("float", [None, input_features]) Y = tf.placeholder("float", [None, output_class]) logits = RNN(X,hidden_nodes0,hidden_nodes,input_features,output_class) loss = tf.losses.mean_squared_error(predictions = logits, labels = Y) global_step = tf.Variable(0) learning_rate = tf.train.exponential_decay( alpha, global_step, num_epochs, 0.99, staircase=True) optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate, epsilon = 1e-10).minimize(loss,global_step=global_step) init = tf.global_variables_initializer() # 训练 with tf.Session() as sess: sess.run(init) N = P.shape[0] for epoch in range(num_epochs): total_batch = int(math.ceil(N / batch_size)) indices = np.arange(N) np.random.shuffle(indices) avg_loss = 0 # 迭代训练，顺便计算训练集loss for i in range(total_batch): rand_index = indices[batch_size*i:batch_size*(i+1)] x = P[rand_index] y = T[rand_index] _, cost = sess.run([optimizer, loss], feed_dict={X: x, Y: y}) avg_loss += cost / total_batch #计算测试集的预测值 test_pred = sess.run(logits, feed_dict={X: Pt}) test_pred = test_pred.reshape(-1, output_class) F2=np.mean(np.square((test_pred-Tt))) return F2 # 对超过边界的变量进行去除 def Bounds(x,lb,ub,dim): temp = x.reshape(-1,dim) I = temp < lb temp[I] = lb[I] J = temp > ub temp[J] = ub[J] return temp def SSA(P,T,Pt,Tt): M=1 pop=5 dim=4 P_percent=0.2 #M 迭代次数 #pop 麻雀种群数量 #dim 寻优维度 #P_percent 麻雀在生产者的比例 pNum = round( pop * P_percent )#pNum是生产者 x = np.zeros((pop, dim)) fit = np.zeros((pop,1)) # 种群初始化 for i in range(pop): for j in range(dim): if j==0: x[i][j] = np.random.randint(1,100)#产生1-100的整数 elif j==1: x[i][j] = (0.01-0.001)*np.random.rand()+0.001# 产生0.001-0.01的小数 elif j==2: x[i][j] = np.random.randint(1,100)#产生1-100的整数 elif j==3: x[i][j] = np.random.randint(1,100)#产生1-100的整数 fit[ i ] = fun( x[ i, : ],P,T,Pt,Tt ) pFit = fit.copy() #最佳适应度矩阵 pX = x.copy() #最佳种群位置 fMin=np.min( fit ) #fMin表示全局最优适应值，生产者能量储备水平取决于对个人适应度值的评估 bestI=np.argmin( fit ) bestX = x[bestI, : ].copy() # bestX表示fMin对应的全局最优位置的变量信息 Convergence_curve=np.zeros((M,)) trace=np.zeros((M,dim)) # 初始化收敛曲线 for t in range(M): sortIndex = np.argsort( pFit.reshape(-1,) ).reshape(-1,) # 对麻雀的适应度值进行排序，并取出下标 fmax=np.max( pFit ) # 取出最大的适应度值 B=np.argmax( pFit ) # 取出最大的适应度值得下标 worse= x[B,:].copy() # 最差适应度 r2=np.random.rand() # 预警值 ## 这一部分为发现者（探索者）的位置更新 if r2<0.8:#%预警值较小，说明没有捕食者出现 for i in range(pNum):#r2小于0.8时发现者改变位置 r1=np.random.rand() x[sortIndex[i],:]=pX[sortIndex[i],:]*np.exp(-i/(r1*M)) x[sortIndex[i],:]=boundary(x[sortIndex[i],:]) temp=fun( x[ sortIndex[ i ], : ],P,T,Pt,Tt ) fit[ sortIndex[ i ] ] = temp# 计算新的适应度值 else:#预警值较大，说明有捕食者出现威胁到了种群的安全，需要去其它地方觅食 for i in range(pNum):#r2大于0.8时发现者改变位置 r1=np.random.rand() x[sortIndex[i],:]=pX[sortIndex[i],:] + np.random.normal()*np.ones((1,dim)) x[sortIndex[i],:]=boundary(x[sortIndex[i],:]) fit[ sortIndex[ i ] ] = fun( x[ sortIndex[ i ], : ],P,T,Pt,Tt )# 计算新的适应度值 bestII=np.argmin( fit ) bestXX = x[ bestII, : ].copy() ##这一部分为加入者（追随者）的位置更新 for i in range(pNum + 1, pop):#剩下的个体变化 A=np.floor(np.random.rand(1, dim)*2)*2-1 if i > pop/2: #这个代表这部分麻雀处于十分饥饿的状态（因为它们的能量很低，也是是适应度值很差），需要到其它地方觅食 x[ sortIndex[i ], : ]=np.random.normal()*np.exp((worse-pX[sortIndex[ i ], : ])/(i**2)) else:#这一部分追随者是围绕最好的发现者周围进行觅食，其间也有可能发生食物的争夺，使其自己变成生产者 x[ sortIndex[ i ], : ]=bestXX+np.abs( pX[ sortIndex[ i ], : ]-bestXX).dot(A.T*(A*A.T)**(-1))*np.ones((1,dim)) x[sortIndex[ i ], : ] = boundary( x[ sortIndex[ i ], : ])#判断边界是否超出 fit[ sortIndex[ i ] ] = fun( x[ sortIndex[ i ], : ],P,T,Pt,Tt )#计算适应度值 #这一部分为意识到危险（注意这里只是意识到了危险，不代表出现了真正的捕食者）的麻雀的位置更新 c=random.sample(range(sortIndex.shape[0]),sortIndex.shape[0])#这个的作用是在种群中随机产生其位置（也就是这部分的麻雀位置一开始是随机的，意识到危险了要进行位置移动， b=sortIndex[np.array(c)[0:round(pop*0.2)]].reshape(-1,) for j in range(b.shape[0]): if pFit[sortIndex[b[j]]]> fMin:#处于种群外围的麻雀的位置改变 x[ sortIndex[ b[j] ], : ]=bestX+np.random.normal(1,dim)*(np.abs( pX[ sortIndex[ b[j] ], : ] -bestX)) else: #处于种群中心的麻雀的位置改变 x[ sortIndex[ b[j] ], : ] =pX[ sortIndex[ b[j] ], : ] + (2*np.random.rand()-1)*(np.abs(pX[ sortIndex[ b[j] ], : ]-worse)) / (pFit[sortIndex[b[j]]]-fmax+1e-50) x[ sortIndex[b[j] ], : ] = boundary( x[ sortIndex[b[j]], : ]) fit[ sortIndex[b[j] ] ] = fun( x[ sortIndex[b[j] ], : ],P,T,Pt,Tt )#计算适应度值 # 这部分是最终的最优解更新 for i in range(pop): if fit[ i ] < pFit[ i ] : pFit[ i ] = fit[ i ].copy() pX[ i, : ] = x[ i, : ].copy() if pFit[i ] < fMin: fMin= pFit[ i ,0].copy() bestX = pX[ i, :].copy() trace[t, :]=bestX print(t+1, fMin, bestX) Convergence_curve[t]=fMin return Convergence_curve, bestX, trace def boundary(pop): # 防止粒子跳出范围 pop[0]=int(pop[0]) pop[2]=int(pop[2]) pop[3]=int(pop[3])#迭代数和节点数都应为整数 if pop[0]>100 or pop[0]<1: pop[0]=np.random.randint(1,PNUMBER) if pop[1]>0.01 or pop[1]<0.001: pop[1] = (0.01-0.001)*np.random.rand()+0.001 if pop[2]>100 or pop[2]<1: pop[2] = np.random.randint(1,100) if pop[3]>100 or pop[3]<1: pop[3] = np.random.randint(1,100) return pop # In[] 加载数据 def split_data(data,n): # 前n时刻的值为输入 # 预测n+1个时刻 in_=[] out_=[] N=data.shape[0]-n for i in range(N): in_.append(data[i:i+n]) out_.append(data[i+n]) in_=np.array(in_).reshape(len(in_),-1) out_=np.array(out_).reshape(len(out_),-1) return in_,out_ def train_one(name,df, num_epochs = 50): # -*- coding: utf-8 -*- # lstm data = df[name].values data = data.reshape(-1,1) data.shape data=data.reshape(-1, 1) in_, out_ = split_data(data, 48) n = range(in_.shape[0]) m = -48#测试集 train_data = in_[n[0:m],] print(train_data.shape) test_data = in_[n[m:],] train_label = out_[n[0:m],] print(train_label.shape) test_label = out_[n[m:],] # 归一化 # ss_X=StandardScaler().fit(train_data) # ss_Y=StandardScaler().fit(train_label) ss_X=MinMaxScaler(feature_range=(0, 1)).fit(train_data) ss_Y=MinMaxScaler(feature_range=(0, 1)).fit(train_label) train_data = ss_X.transform(train_data) train_label = ss_Y.transform(train_label) test_data = ss_X.transform(test_data) print("test_data的shape:",test_data.shape) test_label = ss_Y.transform(test_label) in_num=train_data.shape[1] out_num=train_label.shape[1] trace, best, result = SSA(train_data,train_label,test_data,test_label) ################################################### 取出 num_epochs = int(best[0])#迭代次数 alpha = best[1]# 学习率 hidden_nodes0 = int(best[2])#第一隐含层神经元 hidden_nodes = int(best[3])#第二隐含层神经元 # In[]定义超参数 batch_size = 128# batchsize input_features = in_num output_class = out_num print(in_num) print(out_num) # placeholder X = tf.placeholder("float", [None, input_features]) Y = tf.placeholder("float", [None, output_class]) # In[] 初始化 logits = single_RNN(X, hidden_nodes0, hidden_nodes, input_features,output_class) loss = tf.losses.mean_squared_error(predictions=logits, labels=Y) global_step = tf.Variable(0) learning_rate = tf.train.exponential_decay( alpha, global_step, num_epochs, 0.98, staircase=True) optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate, epsilon=1e-10).minimize(loss, global_step=global_step) init = tf.global_variables_initializer() # In[]训练 train = [] valid=[] with tf.Session() as sess: sess.run(init) N = train_data.shape[0] for epoch in range(num_epochs): total_batch = int(math.ceil(N / batch_size)) indices = np.arange(N) np.random.shuffle(indices) avg_loss = 0 # 迭代训练，计算训练集loss for i in range(total_batch): rand_index = indices[batch_size*i:batch_size*(i+1)] x = train_data[rand_index] y = train_label[rand_index] _, cost = sess.run([optimizer, loss], feed_dict={X: x, Y: y}) avg_loss += cost / total_batch # 计算测试集loss valid_data=test_data.reshape(-1, input_features) valid_y = test_label.reshape(-1, output_class) valid_loss = sess.run(loss, feed_dict={X: valid_data, Y: valid_y}) train.append(avg_loss) valid.append(valid_loss) print('epoch:',epoch,' ,train loss ',avg_loss,' ,valid loss: ',valid_loss) test_data=test_data.reshape(-1,input_features) print(test_data.shape) test_pred = sess.run(logits, feed_dict={X: test_data}) print(test_pred.shape) test_pred = test_pred.reshape(-1, output_class) #结果反归一化 test_label = ss_Y.inverse_transform(test_label) test_pred = ss_Y.inverse_transform(test_pred) return test_label,test_pred,trace